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