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