]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMUONAlignment.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New versions of GDC and CDH raw data headers. Some CDH getters are added
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignment.cxx
CommitLineData
010eb601 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//-----------------------------------------------------------------------------
19/// \class AliMUONAlignmentV5
20/// Alignment class fro the ALICE DiMuon spectrometer
21///
22/// MUON specific alignment class which interface to AliMillepede.
23/// For each track ProcessTrack calculates the local and global derivatives
24/// at each hit and fill the corresponding local equations. Provide methods for
25/// fixing or constraining detection elements for best results.
26///
27/// \author Bruce Becker, Javier Castillo
28//-----------------------------------------------------------------------------
29
30#include "AliMUONAlignment.h"
31#include "AliMUONTrack.h"
32#include "AliMUONRawCluster.h"
33#include "AliMUONTrackParam.h"
34#include "AliMUONHitForRec.h"
35#include "AliMUONGeometryTransformer.h"
36#include "AliMUONGeometryModuleTransformer.h"
37#include "AliMUONGeometryDetElement.h"
38#include "AliMUONGeometryStore.h"
39#include "AliMUONGeometryBuilder.h"
40#include "AliLog.h"
41#include "TSystem.h"
42#include "AliMUONConstants.h"
43
44#include "AliMillepede.h"
45
46ClassImp(AliMUONAlignment)
47 Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
48 Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
49 Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
50 Int_t AliMUONAlignment::fgNParCh = 3;
51 Int_t AliMUONAlignment::fgNCh = 10;
52 Int_t AliMUONAlignment::fgNSt = 5;
53
54AliMUONAlignment::AliMUONAlignment() : TObject()
55{
56 /// Default constructor, setting define alignment parameters
57 fSigma[0] = 1.0e-1;
58 fSigma[1] = 1.0e-2;
59 fPhi = 0;
60 fNGlobal = fgNDetElem*fgNParCh;
61 fNLocal = 4;
62 fNStdDev = 3;
63
64 fDoF[0] = kTRUE; fDoF[1] = kTRUE; fDoF[2] = kTRUE;
65 fAllowVar[0] = 0.05; fAllowVar[1] = 0.05; fAllowVar[2] = 0.001;
66 fStartFac = 16.;
67 fResCut = 100.;
68 fResCutInitial = 100.;
69
70 fBFieldOn = kTRUE;
71// fBFieldOn = kFALSE;
72
73 AliInfo(Form("fAllowVar[0]: %f\t fAllowVar[1]: %f\t fPhi: %f\t fgNDetElem: %i\t fNGlobal: %i\t fNLocal: %i",fAllowVar[0],fAllowVar[1],fPhi,fgNDetElem,fNGlobal,fNLocal));
74
75 fMillepede = new AliMillepede();
76
77 Init(fNGlobal, fNLocal, fNStdDev);
78
79 ResetLocalEquation();
80 AliInfo("Parameters initialized to zero");
81
82}
83
84AliMUONAlignment::~AliMUONAlignment() {
85 /// Destructor
86}
87
88void AliMUONAlignment::Init(Int_t nGlobal, /* number of global paramers */
89 Int_t nLocal, /* number of local parameters */
90 Int_t nStdDev /* std dev cut */ )
91{
92 /// Initialization of AliMillepede. Fix parameters, define constraints ...
93 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
94
95 Bool_t bStOnOff[5] = {kFALSE,kFALSE,kTRUE,kTRUE,kTRUE};
96
97 AllowVariations(bStOnOff);
98
99 // Fix parameters or add constraints here
100 for (Int_t iSt=0; iSt<5; iSt++)
101 if (!bStOnOff[iSt]) FixStation(iSt+1);
102
103 Bool_t bVarXYT[3] = {kFALSE,kTRUE,kFALSE};
104 Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
105 ResetConstraints();
106 AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
107 bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
108 bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
109 AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
110 bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
111 AddConstraints(bStOnOff,bVarXYT);
112
113 // Set iterations
114 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
115}
116
117void AliMUONAlignment::FixStation(Int_t iSt){
118 /// Fix all detection elements of station iSt
119 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
120 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
121 for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){
122 FixParameter(i*fgNParCh+0, 0.0);
123 FixParameter(i*fgNParCh+1, 0.0);
124 FixParameter(i*fgNParCh+2, 0.0);
125 }
126}
127
128void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){
129 /// Set non linear parameter flag selected stations and degrees of freedom
130 for (Int_t i = 0; i < fgNDetElem; i++){
131 Int_t iCh=0;
132 for (iCh=1; iCh<=fgNCh; iCh++){
133 if (i<fgSNDetElemCh[iCh-1]) break;
134 }
135 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
136 if (iSt){
137 if (lVarXYT[0]) { // X constraint
138 SetNonLinear(i*fgNParCh+0);
139 }
140 if (lVarXYT[1]) { // Y constraint
141 SetNonLinear(i*fgNParCh+1);
142 }
143 if (lVarXYT[2]) { // T constraint
144 SetNonLinear(i*fgNParCh+2);
145 }
146 }
147 }
148}
149
150void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT){
151 /// Add constraint equations for selected stations and degrees of freedom
152 for (Int_t i = 0; i < fgNDetElem; i++){
153 Int_t iCh=0;
154 for (iCh=1; iCh<=fgNCh; iCh++){
155 if (i<fgSNDetElemCh[iCh-1]) break;
156 }
157 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
158 if (iSt){
159 if (lVarXYT[0]) { // X constraint
160 fConstraintX[i*fgNParCh+0]=1.0;
161 }
162 if (lVarXYT[1]) { // Y constraint
163 fConstraintY[i*fgNParCh+1]=1.0;
164 }
165 if (lVarXYT[2]) { // T constraint
166 fConstraintP[i*fgNParCh+2]=1.0;
167 }
168 }
169 }
170 if (lVarXYT[0]) { // X constraint
171 AddConstraint(fConstraintX,0.0);
172 }
173 if (lVarXYT[1]) { // Y constraint
174 AddConstraint(fConstraintY,0.0);
175 }
176 if (lVarXYT[2]) { // T constraint
177 AddConstraint(fConstraintP,0.0);
178 }
179}
180
181void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t *lDetTLBR){
182 /// Add constraint equations for selected stations, degrees of freedom detector half
183 for (Int_t i = 0; i < fgNDetElem; i++){
184 Int_t iCh=0;
185 for (iCh=1; iCh<=fgNCh; iCh++){
186 if (i<fgSNDetElemCh[iCh-1]) break;
187 }
188 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
189 if (iSt){
190 if (lVarXYT[0]) { // X constraint
191 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
192 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
193 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
194 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
195 }
196 if (lVarXYT[1]) { // X constraint
197 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
198 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
199 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
200 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
201 }
202 if (lVarXYT[2]) { // X constraint
203 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
204 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
205 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
206 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
207 }
208 }
209 }
210 if (lVarXYT[0]) { // X constraint
211 if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
212 if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
213 if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
214 if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
215 }
216 if (lVarXYT[1]) { // Y constraint
217 if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
218 if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
219 if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
220 if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
221 }
222 if (lVarXYT[2]) { // T constraint
223 if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
224 if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
225 if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
226 if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
227 }
228}
229
230void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar){
231 /// Set constrain equation for top half of spectrometer
232 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
233 if (lCh>=1 && lCh<=4){
234 if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
235 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
236 }
237 }
238 if (lCh>=5 && lCh<=6){
239 if (lDetElemNumber>=0&&lDetElemNumber<=9){
240 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
241 }
242 }
243 if (lCh>=7 && lCh<=10){
244 if (lDetElemNumber>=0&&lDetElemNumber<=13){
245 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
246 }
247 }
248}
249
250void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar){
251 /// Set constrain equation for left half of spectrometer
252 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
253 if (lCh>=1 && lCh<=4){
254 if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
255 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
256 }
257 }
258 if (lCh>=5 && lCh<=6){
259 if (lDetElemNumber>=5&&lDetElemNumber<=13){
260 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
261 }
262 }
263 if (lCh>=7 && lCh<=10){
264 if (lDetElemNumber>=7&&lDetElemNumber<=19){
265 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
266 }
267 }
268}
269
270void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar){
271 /// Set constrain equation for bottom half of spectrometer
272 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
273 if (lCh>=1 && lCh<=4){
274 if (lDetElemNumber==2 && lDetElemNumber==3){ // From track crossings
275 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
276 }
277 }
278 if (lCh>=5 && lCh<=6){
279 if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
280 (lDetElemNumber==0)){
281 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
282 }
283 }
284 if (lCh>=7 && lCh<=10){
285 if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
286 (lDetElemNumber==0)){
287 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
288 }
289 }
290}
291
292void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar){
293 /// Set constrain equation for right half of spectrometer
294 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
295 if (lCh>=1 && lCh<=4){
296 if (lDetElemNumber==0 && lDetElemNumber==3){ // From track crossings
297 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
298 }
299 }
300 if (lCh>=5 && lCh<=6){
301 if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
302 (lDetElemNumber>=14&&lDetElemNumber<=17)){
303 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
304 }
305 }
306 if (lCh>=7 && lCh<=10){
307 if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
308 (lDetElemNumber>=20&&lDetElemNumber<=25)){
309 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
310 }
311 }
312}
313
314void AliMUONAlignment::ResetConstraints(){
315 /// Reset all constraint equations
316 for (Int_t i = 0; i < fgNDetElem; i++){
317 fConstraintX[i*fgNParCh+0]=0.0;
318 fConstraintX[i*fgNParCh+1]=0.0;
319 fConstraintX[i*fgNParCh+2]=0.0;
320 fConstraintY[i*fgNParCh+0]=0.0;
321 fConstraintY[i*fgNParCh+1]=0.0;
322 fConstraintY[i*fgNParCh+2]=0.0;
323 fConstraintP[i*fgNParCh+0]=0.0;
324 fConstraintP[i*fgNParCh+1]=0.0;
325 fConstraintP[i*fgNParCh+2]=0.0;
326 fConstraintXT[i*fgNParCh+0]=0.0;
327 fConstraintXT[i*fgNParCh+1]=0.0;
328 fConstraintXT[i*fgNParCh+2]=0.0;
329 fConstraintYT[i*fgNParCh+0]=0.0;
330 fConstraintYT[i*fgNParCh+1]=0.0;
331 fConstraintYT[i*fgNParCh+2]=0.0;
332 fConstraintPT[i*fgNParCh+0]=0.0;
333 fConstraintPT[i*fgNParCh+1]=0.0;
334 fConstraintPT[i*fgNParCh+2]=0.0;
335 fConstraintXL[i*fgNParCh+0]=0.0;
336 fConstraintXL[i*fgNParCh+1]=0.0;
337 fConstraintXL[i*fgNParCh+2]=0.0;
338 fConstraintYL[i*fgNParCh+0]=0.0;
339 fConstraintYL[i*fgNParCh+1]=0.0;
340 fConstraintYL[i*fgNParCh+2]=0.0;
341 fConstraintPL[i*fgNParCh+0]=0.0;
342 fConstraintPL[i*fgNParCh+1]=0.0;
343 fConstraintPL[i*fgNParCh+2]=0.0;
344 fConstraintXB[i*fgNParCh+0]=0.0;
345 fConstraintXB[i*fgNParCh+1]=0.0;
346 fConstraintXB[i*fgNParCh+2]=0.0;
347 fConstraintYB[i*fgNParCh+0]=0.0;
348 fConstraintYB[i*fgNParCh+1]=0.0;
349 fConstraintYB[i*fgNParCh+2]=0.0;
350 fConstraintPB[i*fgNParCh+0]=0.0;
351 fConstraintPB[i*fgNParCh+1]=0.0;
352 fConstraintPB[i*fgNParCh+2]=0.0;
353 fConstraintXR[i*fgNParCh+0]=0.0;
354 fConstraintXR[i*fgNParCh+1]=0.0;
355 fConstraintXR[i*fgNParCh+2]=0.0;
356 fConstraintYR[i*fgNParCh+0]=0.0;
357 fConstraintYR[i*fgNParCh+1]=0.0;
358 fConstraintYR[i*fgNParCh+2]=0.0;
359 fConstraintPR[i*fgNParCh+0]=0.0;
360 fConstraintPR[i*fgNParCh+1]=0.0;
361 fConstraintPR[i*fgNParCh+2]=0.0;
362 }
363}
364
365void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
366 /// Constrain equation defined by par to value
367 fMillepede->SetGlobalConstraint(par, value);
368 AliInfo("Adding constraint");
369}
370
371void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
372 /// Initialize global parameters with par array
373 fMillepede->SetGlobalParameters(par);
374 AliInfo("Init Global Parameters");
375}
376
377void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
378 /// Parameter iPar is encourage to vary in [-value;value].
379 /// If value == 0, parameter is fixed
380 fMillepede->SetParSigma(iPar, value);
381 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
382}
383
384void AliMUONAlignment::ResetLocalEquation()
385{
386 /// Reset the derivative vectors
387 for(int i=0; i<fNLocal; i++) {
388 fLocalDerivatives[i] = 0.0;
389 }
390 for(int i=0; i<fNGlobal; i++) {
391 fGlobalDerivatives[i] = 0.0;
392 }
393}
394
395void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) {
396 /// Set allowed variation for selected stations based on fDoF and fAllowVar
397 for (Int_t iSt=1; iSt<=5; iSt++) {
398 if (bStOnOff[iSt-1]) {
399 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
400 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
401 for (int i=0; i<fgNParCh; i++) {
402 AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
403 if (fDoF[i]) {
404 for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
405 FixParameter(j*fgNParCh+i, fAllowVar[i]);
406 }
407 }
408 }
409 }
410 }
411}
412
413void AliMUONAlignment::SetNonLinear(Int_t iPar /* set non linear flag */ ) {
414 /// Set nonlinear flag for parameter iPar
415 fMillepede->SetNonLinear(iPar);
416 AliInfo(Form("Parameter %i set to non linear", iPar));
417}
418
419void AliMUONAlignment::LocalEquationX() {
420 /// Define local equation for current track and hit in x coor. measurement
421 // set local derivatives
422 SetLocalDerivative(0, fCosPhi);
423 SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
424 SetLocalDerivative(2, fSinPhi);
425 SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
426
427 // set global derivatives
428 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -1.);
429 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, 0.);
430 if (fBFieldOn){
431 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
432 -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
433 +fCosPhi*(fTrackPos[1]-fDetElemPos[1]));
434 }
435 else {
436 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
437 -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*
438 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
439 +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*
440 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
441 }
442
443 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
444}
445
446void AliMUONAlignment::LocalEquationY() {
447 /// Define local equation for current track and hit in y coor. measurement
448 // set local derivatives
449 SetLocalDerivative(0,-fSinPhi);
450 SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
451 SetLocalDerivative(2, fCosPhi);
452 SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
453
454 // set global derivatives
455 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, 0.);
456 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -1.);
457 if (fBFieldOn){
458 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
459 -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
460 -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
461 }
462 else {
463 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
464 -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*
465 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
466 -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*
467 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
468 }
469
470 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
471}
472
473void AliMUONAlignment::FillRecPointData() {
474 /// Get information of current hit
475 fClustPos[0] = fRecHit->GetNonBendingCoor();
476 fClustPos[1] = fRecHit->GetBendingCoor();
477 fClustPos[2] = fRecHit->GetZ();
478 fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
479 fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
480}
481
482void AliMUONAlignment::FillTrackParamData() {
483 /// Get information of current track at current hit
484 fTrackPos[0] = fTrackParam->GetNonBendingCoor();
485 fTrackPos[1] = fTrackParam->GetBendingCoor();
486 fTrackPos[2] = fTrackParam->GetZ();
487 fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
488 fTrackSlope[1] = fTrackParam->GetBendingSlope();
489 fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
490 fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
491}
492
493void AliMUONAlignment::FillDetElemData() {
494 /// Get information of current detection element
495 Double_t lDetElemLocX = 0.;
496 Double_t lDetElemLocY = 0.;
497 Double_t lDetElemLocZ = 0.;
498 fDetElemId = fRecHit->GetDetElemId();
499 fDetElemNumber = fDetElemId%100;
500 for (int iCh=0; iCh<fDetElemId/100-1; iCh++){
501 fDetElemNumber += fgNDetElemCh[iCh];
502 }
503 fTransform->Local2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
504 fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
505 if (fDetElemId/100 < 5){
506 fSigma[0] = 3.0e-2;
507 fSigma[1] = 3.0e-2;
508 }
509 else {
510 fSigma[0] = 1.0e-1;
511 fSigma[1] = 1.0e-2;
512 }
513}
514
515void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) {
516 /// Process track; Loop over hits and set local equations
517 fTrack = track;
518 // get tclones arrays.
519 fTrackParamAtHit = fTrack->GetTrackParamAtHit();
520 fHitForRecAtHit = fTrack->GetHitForRecAtHit();
521
522 // get size of arrays
523 Int_t nTrackParam = fTrackParamAtHit->GetEntries();
524 Int_t nHitForRec = fHitForRecAtHit->GetEntries();
525 AliInfo(Form("Number of track param entries : %i ", nTrackParam));
526 AliInfo(Form("Number of hit for rec entries : %i ", nHitForRec));
527
528 for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
529 fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
530 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
531 if (!fRecHit || !fTrackParam) continue;
532 // fill local variables for this position --> one measurement
533 FillDetElemData();
534 FillRecPointData();
535 FillTrackParamData();
536// if (fDetElemId<500) continue;
537 fTrackPos0[0] = fTrackPos[0];
538 fTrackPos0[1] = fTrackPos[1];
539 fTrackPos0[2] = fTrackPos[2];
540 fTrackSlope0[0] = fTrackSlope[0];
541 fTrackSlope0[1] = fTrackSlope[1];
542 break;
543 }
544
545 for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
546 // and get new pointers
547 fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
548 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
549 if (!fRecHit || !fTrackParam) continue;
550 // fill local variables for this position --> one measurement
551 FillDetElemData();
552 FillRecPointData();
553 FillTrackParamData();
554// if (fDetElemId<500) continue;
555 AliDebug(1,Form("cluster: %i", iHit));
556 AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
557 AliDebug(1,Form("fDetElemPos[0]: %f\t fDetElemPos[1]: %f\t fDetElemPos[2]: %f\t DetElemID: %i\t ", fDetElemPos[0],fDetElemPos[1],fDetElemPos[2], fDetElemId));
558
559 AliDebug(1,Form("Track Parameter: %i", iHit));
560 AliDebug(1,Form("x: %f\t y: %f\t z: %f\t slopex: %f\t slopey: %f", fTrackPos[0], fTrackPos[1], fTrackPos[2], fTrackSlope[0], fTrackSlope[1]));
561 AliDebug(1,Form("x0: %f\t y0: %f\t z0: %f\t slopex0: %f\t slopey0: %f", fTrackPos0[0], fTrackPos0[1], fTrackPos0[2], fTrackSlope0[0], fTrackSlope0[1]));
562
563 fCosPhi = TMath::Cos(fPhi);
564 fSinPhi = TMath::Sin(fPhi);
565 if (fBFieldOn){
566 fMeas[0] = fTrackPos[0] - fClustPos[0];
567 fMeas[1] = fTrackPos[1] - fClustPos[1];
568 }
569 else {
570 fMeas[0] = - fClustPos[0];
571 fMeas[1] = - fClustPos[1];
572 }
573 AliDebug(1,Form("fMeas[0]: %f\t fMeas[1]: %f\t fSigma[0]: %f\t fSigma[1]: %f", fMeas[0], fMeas[1], fSigma[0], fSigma[1]));
574 // Set local equations
575 LocalEquationX();
576 LocalEquationY();
577 }
578}
579
580void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
581 /// Call local fit for this tracks
582 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
583 if (iRes && !lSingleFit) {
584 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
585 }
586}
587
588void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
589 /// Call global fit; Global parameters are stored in parameters
590 fMillepede->GlobalFit(parameters,errors,pulls);
591 AliInfo("Done fitting global parameters!");
592 for (int iGlob=0; iGlob<fgNDetElem; iGlob++){
593 printf("%d\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+2]);
594 }
595}
596
597Double_t AliMUONAlignment::GetParError(Int_t iPar) {
598 /// Get error of parameter iPar
599 Double_t lErr = fMillepede->GetParError(iPar);
600 return lErr;
601}
602
603void AliMUONAlignment::PrintGlobalParameters() {
604 /// Print global parameters
605 fMillepede->PrintGlobalParameters();
606}
607
608//_________________________________________________________________________
609TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *detElemMisAlignment) const
610{
611 /// Realign given transformation by given misalignment and return the misaligned transformation
612
613 Double_t cartMisAlig[3] = {0,0,0};
614 Double_t angMisAlig[3] = {0,0,0};
615 const Double_t *trans = transform.GetTranslation();
616 TGeoRotation *rot;
617 // check if the rotation we obtain is not NULL
618 if (transform.GetRotation()) {
619 rot = transform.GetRotation();
620 }
621 else {
622 rot = new TGeoRotation("rot");
623 } // default constructor.
624
625 cartMisAlig[0] = -detElemMisAlignment[0];
626 cartMisAlig[1] = -detElemMisAlignment[1];
627 angMisAlig[2] = -detElemMisAlignment[2]*180./TMath::Pi();
628
629 TGeoTranslation newTrans(cartMisAlig[0] + trans[0], cartMisAlig[1] + trans[1], cartMisAlig[2] + trans[2]);
630
631 rot->RotateX(angMisAlig[0]);
632 rot->RotateY(angMisAlig[1]);
633 rot->RotateZ(angMisAlig[2]);
634
635 return TGeoCombiTrans(newTrans, *rot);
636}
637
638//______________________________________________________________________
639AliMUONGeometryTransformer *
640AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
641 double *misAlignments, Bool_t verbose)
642
643{
644 /////////////////////////////////////////////////////////////////////
645 // Takes the internal geometry module transformers, copies them
646 // and gets the Detection Elements from them.
647 // Takes misalignment parameters and applies these
648 // to the local transform of the Detection Element
649 // Obtains the global transform by multiplying the module transformer
650 // transformation with the local transformation
651 // Applies the global transform to a new detection element
652 // Adds the new detection element to a new module transformer
653 // Adds the new module transformer to a new geometry transformer
654 // Returns the new geometry transformer
655
656 Double_t lDetElemMisAlignment[3] = {0.,0.,0.};
657
658 AliMUONGeometryTransformer *newGeometryTransformer =
659 new AliMUONGeometryTransformer(kTRUE);
660 for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++)
661 { // module transformers
662
663 const AliMUONGeometryModuleTransformer *kModuleTransformer =
664 transformer->GetModuleTransformer(iMt, true);
665
666 AliMUONGeometryModuleTransformer *newModuleTransformer =
667 new AliMUONGeometryModuleTransformer(iMt);
668 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
669
670 TGeoCombiTrans moduleTransform =
671 TGeoCombiTrans(*kModuleTransformer->GetTransformation());
672 TGeoCombiTrans *newModuleTransform = new TGeoCombiTrans(moduleTransform);
673 // same module transform as the previous one
674 // no mis align object created
675 newModuleTransformer->SetTransformation(moduleTransform);
676
677 AliMUONGeometryStore *detElements =
678 kModuleTransformer->GetDetElementStore();
679
680 if (verbose)
681 AliInfo(Form
682 ("%i DEs in old GeometryStore %i",
683 detElements->GetNofEntries(), iMt));
684 Int_t iBase = !iMt ? 0 : fgSNDetElemCh[iMt-1];
685 for (Int_t iDe = 0; iDe < detElements->GetNofEntries(); iDe++)
686 { // detection elements.
687 AliMUONGeometryDetElement *detElement =
688 (AliMUONGeometryDetElement *) detElements->GetEntry(iDe);
689 if (!detElement)
690 AliFatal("Detection element not found.");
691
692 /// make a new detection element
693 AliMUONGeometryDetElement *newDetElement =
694 new AliMUONGeometryDetElement(detElement->GetId(),
695 detElement->GetVolumePath());
696 for (int i=0; i<fgNParCh; i++) {
697 AliInfo(Form("iMt %i, iBase %i, iDe %i, iPar %i",iMt, iBase, iDe, (iBase+iDe)*fgNParCh+i));
698 lDetElemMisAlignment[i] =
699 (iMt<fgNCh) ? misAlignments[(iBase+iDe)*fgNParCh+i] : 0.;
700 }
701 // local transformation of this detection element.
702 TGeoCombiTrans localTransform
703 = TGeoCombiTrans(*detElement->GetLocalTransformation());
704 TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
705 newDetElement->SetLocalTransformation(newLocalTransform);
706
707 // global transformation
708 TGeoHMatrix newGlobalTransform =
709 AliMUONGeometryBuilder::Multiply(*newModuleTransform,
710 newLocalTransform);
711 newDetElement->SetGlobalTransformation(newGlobalTransform);
712
713 // add this det element to module
714 newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
715 newDetElement);
716 // Get delta transformation:
717 // Tdelta = Tnew * Told.inverse
718 TGeoHMatrix deltaTransform
719 = AliMUONGeometryBuilder::Multiply(
720 newGlobalTransform,
721 detElement->GetGlobalTransformation()->Inverse());
722
723 // Create mis alignment matrix
724 newGeometryTransformer
725 ->AddMisAlignDetElement(detElement->GetId(), deltaTransform);
726 }
727 if (verbose)
728 AliInfo(Form("Added module transformer %i to the transformer", iMt));
729 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
730 }
731 return newGeometryTransformer;
732}
733