]> 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
ATO-97 - addopt to the new naming covnevention for branches in calibration trees
[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
516043cd 16/* $Id: AliMUONAlignment.cxx 51000 2011-08-08 17:58:17Z ivana $ */
010eb601 17
18//-----------------------------------------------------------------------------
c25af45b 19/// \class AliMUONAlignment
56218672 20/// Alignment class for the ALICE DiMuon spectrometer
010eb601 21///
56218672 22/// MUON specific alignment class which interface to AliMillepede.
010eb601 23/// For each track ProcessTrack calculates the local and global derivatives
96ebe67e 24/// at each cluster and fill the corresponding local equations. Provide methods
56218672 25/// for fixing or constraining detection elements for best results.
010eb601 26///
27/// \author Bruce Becker, Javier Castillo
28//-----------------------------------------------------------------------------
29
30#include "AliMUONAlignment.h"
31#include "AliMUONTrack.h"
010eb601 32#include "AliMUONTrackParam.h"
96ebe67e 33#include "AliMUONVCluster.h"
010eb601 34#include "AliMUONGeometryTransformer.h"
35#include "AliMUONGeometryModuleTransformer.h"
36#include "AliMUONGeometryDetElement.h"
010eb601 37#include "AliMUONGeometryBuilder.h"
516043cd 38#include "AliMillePede2.h"
010eb601 39
8c95578a 40#include "AliMpExMap.h"
630711ed 41#include "AliMpExMapIterator.h"
8c95578a 42
4818a9b7 43#include "AliAlignObjMatrix.h"
8c95578a 44#include "AliLog.h"
45
516043cd 46#include <TMath.h>
47#include <TMatrixDSym.h>
48#include <TClonesArray.h>
49#include <TGraphErrors.h>
8c95578a 50
78649106 51/// \cond CLASSIMP
010eb601 52ClassImp(AliMUONAlignment)
78649106 53/// \endcond
54
56218672 55//_____________________________________________________________________
56// static variables
516043cd 57const Int_t AliMUONAlignment::fgNDetElemCh[AliMUONAlignment::fgNCh] = { 4, 4, 4, 4, 18, 18, 26, 26, 26, 26 };
58const Int_t AliMUONAlignment::fgSNDetElemCh[AliMUONAlignment::fgNCh+1] = { 0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156 };
59
60//_____________________________________________________________________
61/// self initialized array, used for adding constraints
62class Array
63{
64
65 public:
66
67 /// contructor
68 Array( void )
69 {
70 for( Int_t i=0; i < AliMUONAlignment::fNGlobal; ++i )
71 { values[i] = 0; }
72 }
73
74 /// array
75 Double_t values[AliMUONAlignment::fNGlobal];
76
77 private:
78
79 /// Not implemented
80 Array(const Array& );
81
82 /// Not implemented
83 Array& operator = (const Array& );
84
85};
56218672 86
87//_____________________________________________________________________
88AliMUONAlignment::AliMUONAlignment()
8395bff1 89 : TObject(),
516043cd 90 fInitialized( kFALSE ),
91 fRunNumber( 0 ),
92 fBFieldOn( kTRUE ),
93 fStartFac( 256 ),
94 fResCutInitial( 100 ),
95 fResCut( 100 ),
96 fMillepede( 0L ),
97 fCluster( 0L ),
98 fNStdDev( 3 ),
99 fDetElemNumber( 0 ),
56218672 100 fTrackRecord(),
516043cd 101 fTransform( 0 ),
102 fGeoCombiTransInverse()
010eb601 103{
516043cd 104 /// constructor
56218672 105 fSigma[0] = 1.5e-1;
010eb601 106 fSigma[1] = 1.0e-2;
010eb601 107
516043cd 108 // default allowed variations
109 fAllowVar[0] = 0.5; // x
110 fAllowVar[1] = 0.5; // y
111 fAllowVar[2] = 0.01; // phi_z
112 fAllowVar[3] = 5; // z
64e2c144 113
516043cd 114 // initialize millepede
115 fMillepede = new AliMillePede2();
56218672 116
516043cd 117 // initialize degrees of freedom
118 // by default all parameters are free
119 for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
120 { fGlobalParameterStatus[iPar] = kFreeParId; }
010eb601 121
516043cd 122 // initialize local equations
123 for(int i=0; i<fNLocal; ++i )
124 { fLocalDerivatives[i] = 0.0; }
010eb601 125
516043cd 126 for(int i=0; i<fNGlobal; ++i )
127 { fGlobalDerivatives[i] = 0.0; }
010eb601 128
129}
130
d9426e3f 131//_____________________________________________________________________
516043cd 132AliMUONAlignment::~AliMUONAlignment()
d9426e3f 133{
516043cd 134 /// destructor
d9426e3f 135}
136
56218672 137//_____________________________________________________________________
516043cd 138void AliMUONAlignment::Init( void )
010eb601 139{
010eb601 140
516043cd 141 /// initialize
142 /**
143 initialize millipede
144 must be called after necessary detectors have been fixed,
145 but before constrains are added and before global parameters initial value are set
146 */
147 if( fInitialized )
148 { AliFatal( "Millepede already initialized" ); }
149
150 // assign proper groupID to free parameters
151 Int_t nGlobal = 0;
152 for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
153 {
010eb601 154
516043cd 155 if( fGlobalParameterStatus[iPar] == kFixedParId )
156 {
157 // fixed parameters are left unchanged
158 continue;
010eb601 159
516043cd 160 } else if( fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId ) {
56218672 161
516043cd 162 // free parameters or first element of group are assigned a new group id
163 fGlobalParameterStatus[iPar] = nGlobal++;
164 continue;
010eb601 165
516043cd 166 } else if( fGlobalParameterStatus[iPar] < kGroupBaseId ) {
3b2890be 167
516043cd 168 // get detector element id from status, get chamber parameter id
169 const Int_t iDeBase( kGroupBaseId - 1 - fGlobalParameterStatus[iPar] );
170 const Int_t iParBase = iPar%fgNParCh;
56218672 171
516043cd 172 // check
173 if( iDeBase < 0 || iDeBase >= iPar/fgNParCh )
174 { AliFatal( Form( "Group for parameter index %i has wrong base detector element: %i", iPar, iDeBase ) ); }
3b2890be 175
516043cd 176 // assign identical group id to current
177 fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase*fgNParCh + iParBase];
178 AliInfo( Form( "Parameter %i grouped to detector %i (%s)", iPar, iDeBase, GetParameterMaskString( 1<<iParBase ).Data() ) );
3b2890be 179
516043cd 180 } else AliFatal( Form( "Unrecognized parameter status for index %i: %i", iPar, fGlobalParameterStatus[iPar] ) );
181
182 }
183
184 AliInfo( Form( "Free Parameters: %i out of %i", nGlobal, fNGlobal ) );
185
186 // initialize millepede
187 fMillepede->InitMille( fNGlobal, fNLocal, fNStdDev, fResCut, fResCutInitial, fGlobalParameterStatus );
188 fInitialized = kTRUE;
189
190 // some debug output
191 for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
192 { AliInfo( Form( "fAllowVar[%i]= %f", iPar, fAllowVar[iPar] ) ); }
193
194 // set allowed variations for all parameters
195 for( Int_t iDet = 0; iDet < fgNDetElem; ++iDet )
196 {
197 for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
198 { fMillepede->SetParSigma( iDet*fgNParCh + iPar, fAllowVar[iPar] ); }
199 }
56218672 200
010eb601 201 // Set iterations
56218672 202 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
516043cd 203
010eb601 204}
205
516043cd 206//_____________________________________________________
207AliMillePedeRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight )
56218672 208{
516043cd 209 /// process track for alignment minimization
210 /**
211 returns the alignment records for this track.
212 They can be stored in some output for later reprocessing.
213 */
214
215 // reset track records
216 fTrackRecord.Reset();
217 if( fMillepede->GetRecord() ) fMillepede->GetRecord()->Reset();
218
219 // get number of track parameters
220 Int_t nTrackParam = track->GetTrackParamAtCluster()->GetEntries();
221
222 Bool_t first( kTRUE );
223 for( Int_t iTrackParam = 0; iTrackParam < nTrackParam; ++iTrackParam )
56218672 224 {
516043cd 225
226 // get new pointers
227 AliMUONTrackParam* trackParam( (AliMUONTrackParam *) track->GetTrackParamAtCluster()->At(iTrackParam) );
228 if( !trackParam ) continue;
229
230 AliMUONVCluster* cluster = trackParam->GetClusterPtr();
231 if( !cluster ) continue;
232
233 // fill local variables for this position --> one measurement
234 FillDetElemData( cluster );
235 FillRecPointData( cluster );
516043cd 236 FillTrackParamData( trackParam );
237
238 if( first )
239 {
240
241 // for first valid cluster, save track position as "starting" values
242 first = kFALSE;
243
244 fTrackPos0[0] = fTrackPos[0];
245 fTrackPos0[1] = fTrackPos[1];
246 fTrackPos0[2] = fTrackPos[2];
247 fTrackSlope0[0] = fTrackSlope[0];
248 fTrackSlope0[1] = fTrackSlope[1];
249
250 }
251
252 // 'inverse' (GlobalToLocal) rotation matrix
253 const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
254
255 // calculate measurements
256 if( fBFieldOn )
257 {
258
259 // use residuals (cluster - track) for measurement
260 fMeas[0] = r[0]*(fClustPos[0] - fTrackPos[0]) + r[1]*(fClustPos[1] - fTrackPos[1]);
261 fMeas[1] = r[3]*(fClustPos[0] - fTrackPos[0]) + r[4]*(fClustPos[1] - fTrackPos[1]);
262
263 } else {
264
265 // use cluster position for measurement
266 fMeas[0] = ( r[0]*fClustPos[0] + r[1]*fClustPos[1] );
267 fMeas[1] = ( r[3]*fClustPos[0] + r[4]*fClustPos[1] );
268
269 }
270
271 // Set local equations
272 LocalEquationX();
273 LocalEquationY();
274
275 }
276
277 // copy track record
278 fMillepede->SetRecordRun(fRunNumber);
279 fMillepede->SetRecordWeight(weight);
280 fTrackRecord = *fMillepede->GetRecord();
281
282 // save record data
283 if( doAlignment )
284 {
285 fMillepede->SaveRecordData();
010eb601 286 }
56218672 287
516043cd 288 // return record
289 return &fTrackRecord;
290
291}
292
293//______________________________________________________________________________
294void AliMUONAlignment::ProcessTrack( AliMillePedeRecord* trackRecord )
295{
296 /// process track record
297 if( !trackRecord ) return;
298
299 // make sure record storage is initialized
300 if( !fMillepede->GetRecord() ) fMillepede->InitDataRecStorage();
301
302 // copy content
303 *fMillepede->GetRecord() = *trackRecord;
304
305 // save record
306 fMillepede->SaveRecordData();
307
308 return;
309
010eb601 310}
311
56218672 312//_____________________________________________________________________
516043cd 313void AliMUONAlignment::FixAll( UInt_t mask )
56218672 314{
516043cd 315 /// fix parameters matching mask, for all chambers
316 AliInfo( Form( "Fixing %s for all detector elements", GetParameterMaskString( mask ).Data() ) );
317
318 // fix all stations
319 for( Int_t i = 0; i < fgNDetElem; ++i )
56218672 320 {
516043cd 321 if( mask & ParX ) FixParameter(i, 0);
322 if( mask & ParY ) FixParameter(i, 1);
323 if( mask & ParTZ ) FixParameter(i, 2);
324 if( mask & ParZ ) FixParameter(i, 3);
723b0b5b 325 }
516043cd 326
723b0b5b 327}
328
56218672 329//_____________________________________________________________________
516043cd 330void AliMUONAlignment::FixChamber( Int_t iCh, UInt_t mask )
56218672 331{
516043cd 332 /// fix parameters matching mask, for all detector elements in a given chamber, counting from 1
56218672 333
516043cd 334 // check boundaries
335 if( iCh < 1 || iCh > 10 )
336 { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
56218672 337
516043cd 338 // get first and last element
339 const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
340 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
341 for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
56218672 342 {
56218672 343
516043cd 344 AliInfo( Form( "Fixing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
56218672 345
516043cd 346 if( mask & ParX ) FixParameter(i, 0);
347 if( mask & ParY ) FixParameter(i, 1);
348 if( mask & ParTZ ) FixParameter(i, 2);
349 if( mask & ParZ ) FixParameter(i, 3);
56218672 350
d8ad38e5 351 }
516043cd 352}
353
354//_____________________________________________________________________
355void AliMUONAlignment::FixDetElem( Int_t iDetElemId, UInt_t mask )
356{
357 /// fix parameters matching mask, for a given detector element, counting from 0
358 const Int_t iDet( GetDetElemNumber( iDetElemId ) );
359 if ( mask & ParX ) FixParameter(iDet, 0);
360 if ( mask & ParY ) FixParameter(iDet, 1);
361 if ( mask & ParTZ ) FixParameter(iDet, 2);
362 if ( mask & ParZ ) FixParameter(iDet, 3);
56218672 363
d8ad38e5 364}
365
56218672 366//_____________________________________________________________________
516043cd 367void AliMUONAlignment::FixHalfSpectrometer( const Bool_t *lChOnOff, UInt_t sidesMask, UInt_t mask )
56218672 368{
369
516043cd 370 /// Fix parameters matching mask for all detectors in selected chambers and selected sides of the spectrometer
371 for( Int_t i = 0; i < fgNDetElem; ++i )
56218672 372 {
373
516043cd 374 // get chamber matching detector
375 const Int_t iCh( GetChamberId(i) );
376 if( !lChOnOff[iCh-1] ) continue;
56218672 377
516043cd 378 // get detector element in chamber
379 Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
380
381 // skip detector if its side is off
382 // stations 1 and 2
383 if( iCh>=1 && iCh<=4 )
56218672 384 {
516043cd 385 if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
386 if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
387 if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
388 if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
389 }
56218672 390
516043cd 391 // station 3
392 if (iCh>=5 && iCh<=6)
393 {
394 if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
395 if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
396 if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
397 if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
398 }
56218672 399
516043cd 400 // stations 4 and 5
401 if (iCh>=7 && iCh<=10)
402 {
403 if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
404 if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
405 if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
406 if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
407 }
56218672 408
516043cd 409 // detector is accepted, fix it
410 FixDetElem( i, mask );
56218672 411
516043cd 412 }
56218672 413
516043cd 414}
56218672 415
516043cd 416//______________________________________________________________________
417void AliMUONAlignment::FixParameter( Int_t iPar )
418{
56218672 419
516043cd 420 /// fix a given parameter, counting from 0
421 if( fInitialized )
422 { AliFatal( "Millepede already initialized" ); }
56218672 423
516043cd 424 fGlobalParameterStatus[iPar] = kFixedParId;
56218672 425
516043cd 426}
56218672 427
56218672 428
516043cd 429//_____________________________________________________________________
430void AliMUONAlignment::ReleaseChamber( Int_t iCh, UInt_t mask )
431{
432 /// release parameters matching mask, for all detector elements in a given chamber, counting from 1
56218672 433
516043cd 434 // check boundaries
435 if( iCh < 1 || iCh > 10 )
436 { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
56218672 437
516043cd 438 // get first and last element
439 const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
440 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
441 for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
442 {
56218672 443
516043cd 444 AliInfo( Form( "Releasing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
445
446 if( mask & ParX ) ReleaseParameter(i, 0);
447 if( mask & ParY ) ReleaseParameter(i, 1);
448 if( mask & ParTZ ) ReleaseParameter(i, 2);
449 if( mask & ParZ ) ReleaseParameter(i, 3);
56218672 450
3b2890be 451 }
516043cd 452}
453
454//_____________________________________________________________________
455void AliMUONAlignment::ReleaseDetElem( Int_t iDetElemId, UInt_t mask )
456{
457 /// release parameters matching mask, for a given detector element, counting from 0
458 const Int_t iDet( GetDetElemNumber( iDetElemId ) );
459 if ( mask & ParX ) ReleaseParameter(iDet, 0);
460 if ( mask & ParY ) ReleaseParameter(iDet, 1);
461 if ( mask & ParTZ ) ReleaseParameter(iDet, 2);
462 if ( mask & ParZ ) ReleaseParameter(iDet, 3);
56218672 463
3b2890be 464}
465
56218672 466//______________________________________________________________________
516043cd 467void AliMUONAlignment::ReleaseParameter( Int_t iPar )
56218672 468{
469
516043cd 470 /// release a given parameter, counting from 0
471 if( fInitialized )
472 { AliFatal( "Millepede already initialized" ); }
56218672 473
516043cd 474 fGlobalParameterStatus[iPar] = kFreeParId;
56218672 475
516043cd 476}
56218672 477
516043cd 478//_____________________________________________________________________
479void AliMUONAlignment::GroupChamber( Int_t iCh, UInt_t mask )
480{
481 /// group parameters matching mask for all detector elements in a given chamber, counting from 1
482 if( iCh < 1 || iCh > 10 )
483 { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
56218672 484
516043cd 485 const Int_t detElemMin = 100*iCh;
486 const Int_t detElemMax = 100*iCh + fgNDetElemCh[iCh]-1;
487 GroupDetElems( detElemMin, detElemMax, mask );
56218672 488
516043cd 489}
56218672 490
516043cd 491//_____________________________________________________________________
492void AliMUONAlignment::GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t mask )
493{
494 /// group parameters matching mask for all detector elements between min and max
495 // check number of detector elements
496 const Int_t nDetElem = detElemMax - detElemMin + 1;
497 if( nDetElem<2 )
498 { AliFatal( Form( "Requested group of DEs %d-%d contains less than 2 DE's", detElemMin, detElemMax ) ); }
56218672 499
516043cd 500 // create list
501 Int_t* detElemList = new int[nDetElem];
502 for( Int_t i = 0; i < nDetElem; ++i )
503 { detElemList[i] = detElemMin+i; }
504
505 // group
506 GroupDetElems( detElemList, nDetElem, mask );
507 delete[] detElemList;
508
509}
510
511//_____________________________________________________________________
512void AliMUONAlignment::GroupDetElems( Int_t* detElemList, Int_t nDetElem, UInt_t mask )
513{
514 /// group parameters matching mask for all detector elements in list
515 if( fInitialized )
516 { AliFatal( "Millepede already initialized" ); }
56218672 517
516043cd 518 const Int_t iDeBase( GetDetElemNumber( detElemList[0] ) );
519 for( Int_t i = 0; i < nDetElem; ++i )
520 {
521 const Int_t iDeCurrent( GetDetElemNumber( detElemList[i] ) );
522 if( mask & ParX ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 0] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
523 if( mask & ParY ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 1] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
524 if( mask & ParTZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 2] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
525 if( mask & ParZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 3] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
526
527 if( i== 0 ) AliInfo( Form( "Creating new group for detector %i and variable %s", detElemList[i], GetParameterMaskString( mask ).Data() ) );
528 else AliInfo( Form( "Adding detector element %i to current group", detElemList[i] ) );
010eb601 529 }
56218672 530
010eb601 531}
532
56218672 533//______________________________________________________________________
516043cd 534void AliMUONAlignment::SetChamberNonLinear( Int_t iCh, UInt_t mask )
56218672 535{
516043cd 536 /// Set parameters matching mask as non linear, for all detector elements in a given chamber, counting from 1
537 const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
538 const Int_t iDetElemLast = fgSNDetElemCh[iCh];
539 for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
56218672 540 {
541
516043cd 542 if( mask & ParX ) SetParameterNonLinear(i, 0);
543 if( mask & ParY ) SetParameterNonLinear(i, 1);
544 if( mask & ParTZ ) SetParameterNonLinear(i, 2);
545 if( mask & ParZ ) SetParameterNonLinear(i, 3);
56218672 546
516043cd 547 }
56218672 548
516043cd 549}
56218672 550
516043cd 551//_____________________________________________________________________
552void AliMUONAlignment::SetDetElemNonLinear( Int_t iDetElemId, UInt_t mask )
553{
554 /// Set parameters matching mask as non linear, for a given detector element, counting from 0
555 const Int_t iDet( GetDetElemNumber( iDetElemId ) );
556 if ( mask & ParX ) SetParameterNonLinear(iDet, 0);
557 if ( mask & ParY ) SetParameterNonLinear(iDet, 1);
558 if ( mask & ParTZ ) SetParameterNonLinear(iDet, 2);
559 if ( mask & ParZ ) SetParameterNonLinear(iDet, 3);
56218672 560
516043cd 561}
56218672 562
516043cd 563//______________________________________________________________________
564void AliMUONAlignment::SetParameterNonLinear( Int_t iPar )
565{
566 /// Set nonlinear flag for parameter iPar
567 if( !fInitialized )
568 { AliFatal( "Millepede not initialized" ); }
56218672 569
516043cd 570 fMillepede->SetNonLinear( iPar );
571 AliInfo( Form( "Parameter %i set to non linear", iPar ) );
572}
573
574//______________________________________________________________________
575void AliMUONAlignment::AddConstraints( const Bool_t *lChOnOff, UInt_t mask )
576{
577 /// Add constraint equations for selected chambers and degrees of freedom
578
579 Array fConstraintX;
580 Array fConstraintY;
581 Array fConstraintTZ;
582 Array fConstraintZ;
56218672 583
516043cd 584 for( Int_t i = 0; i < fgNDetElem; ++i )
56218672 585 {
516043cd 586
587 // get chamber matching detector
588 const Int_t iCh( GetChamberId(i) );
589 if (lChOnOff[iCh-1])
590 {
591
592 if( mask & ParX ) fConstraintX.values[i*fgNParCh+0]=1.0;
593 if( mask & ParY ) fConstraintY.values[i*fgNParCh+1]=1.0;
594 if( mask & ParTZ ) fConstraintTZ.values[i*fgNParCh+2]=1.0;
595 if( mask & ParZ ) fConstraintTZ.values[i*fgNParCh+3]=1.0;
596
597 }
010eb601 598 }
56218672 599
516043cd 600 if( mask & ParX ) AddConstraint(fConstraintX.values,0.0);
601 if( mask & ParY ) AddConstraint(fConstraintY.values,0.0);
602 if( mask & ParTZ ) AddConstraint(fConstraintTZ.values,0.0);
603 if( mask & ParZ ) AddConstraint(fConstraintZ.values,0.0);
604
010eb601 605}
606
56218672 607//______________________________________________________________________
516043cd 608void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask )
56218672 609{
516043cd 610 /*
611 questions:
612 - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
613 - why is weight ignored for ConstrainT and ConstrainB
614 - why is there no constrain on z
615 */
616
56218672 617 /// Add constraint equations for selected chambers, degrees of freedom and detector half
3b2890be 618 Double_t lMeanY = 0.;
619 Double_t lSigmaY = 0.;
620 Double_t lMeanZ = 0.;
621 Double_t lSigmaZ = 0.;
622 Int_t lNDetElem = 0;
516043cd 623
624 for( Int_t i = 0; i < fgNDetElem; ++i )
56218672 625 {
626
516043cd 627 // get chamber matching detector
628 const Int_t iCh( GetChamberId(i) );
629
630 // skip detector if chamber is off
631 if( lChOnOff[iCh-1] ) continue;
632
633 // get detector element id from detector element number
634 const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
635 const Int_t lDetElemId = iCh*100+lDetElemNumber;
636
637 // skip detector if its side is off
638 // stations 1 and 2
639 if( iCh>=1 && iCh<=4 )
640 {
641 if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
642 if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
643 if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
644 if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
010eb601 645 }
516043cd 646
647 // station 3
648 if (iCh>=5 && iCh<=6)
649 {
650 if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
651 if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
652 if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
653 if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
654 }
655
656 // stations 4 and 5
657 if (iCh>=7 && iCh<=10)
658 {
659 if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
660 if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
661 if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
662 if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
663 }
664
665 // get global x, y and z position
666 Double_t lDetElemGloX = 0.;
667 Double_t lDetElemGloY = 0.;
668 Double_t lDetElemGloZ = 0.;
669 fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
670
671 // increment mean Y, mean Z, sigmas and number of accepted detectors
56218672 672 lMeanY += lDetElemGloY;
673 lSigmaY += lDetElemGloY*lDetElemGloY;
674 lMeanZ += lDetElemGloZ;
675 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
676 lNDetElem++;
516043cd 677
56218672 678 }
516043cd 679
680 // calculate mean values
681 lMeanY /= lNDetElem;
682 lSigmaY /= lNDetElem;
683 lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
684 lMeanZ /= lNDetElem;
685 lSigmaZ /= lNDetElem;
686 lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
687 AliInfo( Form( "Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ ) );
688
689 // create all possible arrays
690 Array fConstraintX[4]; //Array for constraint equation X
691 Array fConstraintY[4]; //Array for constraint equation Y
692 Array fConstraintP[4]; //Array for constraint equation P
693 Array fConstraintXZ[4]; //Array for constraint equation X vs Z
694 Array fConstraintYZ[4]; //Array for constraint equation Y vs Z
695 Array fConstraintPZ[4]; //Array for constraint equation P vs Z
696
697 // do we really need these ?
698 Array fConstraintXY[4]; //Array for constraint equation X vs Y
699 Array fConstraintYY[4]; //Array for constraint equation Y vs Y
700 Array fConstraintPY[4]; //Array for constraint equation P vs Y
701
702 // fill Bool_t sides array based on masks, for convenience
703 Bool_t lDetTLBR[4];
704 lDetTLBR[0] = sidesMask & SideTop;
705 lDetTLBR[1] = sidesMask & SideLeft;
706 lDetTLBR[2] = sidesMask & SideBottom;
707 lDetTLBR[3] = sidesMask & SideRight;
708
709 for( Int_t i = 0; i < fgNDetElem; ++i )
710 {
711
712 // get chamber matching detector
713 const Int_t iCh( GetChamberId(i) );
714
715 // skip detector if chamber is off
716 if( !lChOnOff[iCh-1] ) continue;
717
718 // get detector element id from detector element number
719 const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
720 const Int_t lDetElemId = iCh*100+lDetElemNumber;
721
722 // get global x, y and z position
723 Double_t lDetElemGloX = 0.;
724 Double_t lDetElemGloY = 0.;
725 Double_t lDetElemGloZ = 0.;
726 fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
727
728 // loop over sides
729 for( Int_t iSide = 0; iSide < 4; iSide++ )
730 {
731
732 // skip if side is not selected
733 if( !lDetTLBR[iSide] ) continue;
734
735 // skip detector if it is not in the selected side
736 // stations 1 and 2
737 if( iCh>=1 && iCh<=4 )
738 {
739 if( lDetElemNumber == 0 && !(iSide == 0 || iSide == 3) ) continue; // top-right
740 if( lDetElemNumber == 1 && !(iSide == 0 || iSide == 1) ) continue; // top-left
741 if( lDetElemNumber == 2 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
742 if( lDetElemNumber == 3 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
3b2890be 743 }
516043cd 744
745 // station 3
746 if (iCh>=5 && iCh<=6)
747 {
748 if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3) ) continue; // top-right
749 if( lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1) ) continue; // top-left
750 if( lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
751 if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
3b2890be 752 }
516043cd 753
754 // stations 4 and 5
755 if (iCh>=7 && iCh<=10)
756 {
757 if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3) ) continue; // top-right
758 if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1) ) continue; // top-left
759 if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
760 if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
3b2890be 761 }
010eb601 762
516043cd 763 // constrain x
764 if( lVarXYT[0] ) fConstraintX[iSide].values[i*fgNParCh+0] = 1;
765
766 // constrain y
767 if( lVarXYT[1] ) fConstraintY[iSide].values[i*fgNParCh+1] = 1;
768
769 // constrain phi (rotation around z)
770 if( lVarXYT[2] ) fConstraintP[iSide].values[i*fgNParCh+2] = 1;
771
772 // x-z shearing
773 if( lVarXYT[3] ) fConstraintXZ[iSide].values[i*fgNParCh+0] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
774
775 // y-z shearing
776 if( lVarXYT[4] ) fConstraintYZ[iSide].values[i*fgNParCh+1] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
777
778 // phi-z shearing
779 if( lVarXYT[5] ) fConstraintPZ[iSide].values[i*fgNParCh+2] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
780
781 // x-y shearing
782 if( lVarXYT[6] ) fConstraintXY[iSide].values[i*fgNParCh+0] = (lDetElemGloY-lMeanY)/lSigmaY;
783
784 // y-y shearing
785 if( lVarXYT[7] ) fConstraintYY[iSide].values[i*fgNParCh+1] = (lDetElemGloY-lMeanY)/lSigmaY;
786
787 // phi-y shearing
788 if( lVarXYT[8] ) fConstraintPY[iSide].values[i*fgNParCh+2] = (lDetElemGloY-lMeanY)/lSigmaY;
789
010eb601 790 }
516043cd 791
010eb601 792 }
516043cd 793
794 // pass constraints to millepede
795 for( Int_t iSide = 0; iSide < 4; iSide++ )
796 {
797 // skip if side is not selected
798 if( !lDetTLBR[iSide] ) continue;
799
800 if( lVarXYT[0] ) AddConstraint(fConstraintX[iSide].values,0.0);
801 if( lVarXYT[1] ) AddConstraint(fConstraintY[iSide].values,0.0);
802 if( lVarXYT[2] ) AddConstraint(fConstraintP[iSide].values,0.0);
803 if( lVarXYT[3] ) AddConstraint(fConstraintXZ[iSide].values,0.0);
804 if( lVarXYT[4] ) AddConstraint(fConstraintYZ[iSide].values,0.0);
805 if( lVarXYT[5] ) AddConstraint(fConstraintPZ[iSide].values,0.0);
806 if( lVarXYT[6] ) AddConstraint(fConstraintXY[iSide].values,0.0);
807 if( lVarXYT[7] ) AddConstraint(fConstraintYY[iSide].values,0.0);
808 if( lVarXYT[8] ) AddConstraint(fConstraintPY[iSide].values,0.0);
010eb601 809 }
516043cd 810
010eb601 811}
812
56218672 813//______________________________________________________________________
516043cd 814void AliMUONAlignment::InitGlobalParameters(Double_t *par)
815{
816 /// Initialize global parameters with par array
817 if( !fInitialized )
818 { AliFatal( "Millepede is not initialized" ); }
819
820 fMillepede->SetGlobalParameters(par);
010eb601 821}
822
56218672 823//______________________________________________________________________
516043cd 824void AliMUONAlignment::SetAllowedVariation( Int_t iPar, Double_t value )
825{
826 /// "Encouraged" variation for degrees of freedom
827 // check initialization
828 if( fInitialized )
829 { AliFatal( "Millepede already initialized" ); }
830
831 // check initialization
832 if( !(iPar >= 0 && iPar < fgNParCh ) )
833 { AliFatal( Form( "Invalid index: %i", iPar ) ); }
834
835 fAllowVar[iPar] = value;
010eb601 836}
837
56218672 838//______________________________________________________________________
516043cd 839void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
840{
841
842 /// Set expected measurement resolution
843 fSigma[0] = sigmaX;
844 fSigma[1] = sigmaY;
845
846 // print
847 for( Int_t i=0; i<2; ++i )
848 { AliInfo( Form( "fSigma[%i]=%f", i, fSigma[i] ) ); }
849
010eb601 850}
851
516043cd 852//_____________________________________________________
853void AliMUONAlignment::GlobalFit( Double_t *parameters, Double_t *errors, Double_t *pulls )
854{
855
856 /// Call global fit; Global parameters are stored in parameters
857 fMillepede->GlobalFit( parameters, errors, pulls );
858
859 AliInfo( "Done fitting global parameters" );
860 for( int iDet=0; iDet<fgNDetElem; ++iDet )
861 {
862 AliInfo( Form( "%d\t %f\t %f\t %f\t %f",
863 iDet,
864 parameters[iDet*fgNParCh+0],parameters[iDet*fgNParCh+1],
865 parameters[iDet*fgNParCh+3],parameters[iDet*fgNParCh+2]
866 ) );
010eb601 867 }
516043cd 868
010eb601 869}
870
516043cd 871//_____________________________________________________
872void AliMUONAlignment::PrintGlobalParameters() const
873{ fMillepede->PrintGlobalParameters(); }
874
875//_____________________________________________________
876Double_t AliMUONAlignment::GetParError(Int_t iPar) const
877{ return fMillepede->GetParError(iPar); }
878
879//______________________________________________________________________
880AliMUONGeometryTransformer* AliMUONAlignment::ReAlign(
881 const AliMUONGeometryTransformer * transformer,
882 const double *misAlignments, Bool_t )
883{
884
885 /// Returns a new AliMUONGeometryTransformer with the found misalignments
886 /// applied.
887
888 // Takes the internal geometry module transformers, copies them
889 // and gets the Detection Elements from them.
890 // Takes misalignment parameters and applies these
891 // to the local transform of the Detection Element
892 // Obtains the global transform by multiplying the module transformer
893 // transformation with the local transformation
894 // Applies the global transform to a new detection element
895 // Adds the new detection element to a new module transformer
896 // Adds the new module transformer to a new geometry transformer
897 // Returns the new geometry transformer
898
899 Double_t lModuleMisAlignment[fgNParCh] = {0};
900 Double_t lDetElemMisAlignment[fgNParCh] = {0};
901 const TClonesArray* oldMisAlignArray( transformer->GetMisAlignmentData() );
010eb601 902
516043cd 903 AliMUONGeometryTransformer *newGeometryTransformer = new AliMUONGeometryTransformer();
904 for( Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt )
905 {
56218672 906
516043cd 907 // module transformers
908 const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
010eb601 909
516043cd 910 AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
911 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
010eb601 912
516043cd 913 // get transformation
914 TGeoHMatrix deltaModuleTransform( DeltaTransform( lModuleMisAlignment ) );
56218672 915
516043cd 916 // update module
917 TGeoHMatrix moduleTransform( *kModuleTransformer->GetTransformation() );
918 TGeoHMatrix newModuleTransform( AliMUONGeometryBuilder::Multiply( deltaModuleTransform, moduleTransform ) );
919 newModuleTransformer->SetTransformation(newModuleTransform);
920
921 // Get matching old alignment and update current matrix accordingly
922 if( oldMisAlignArray )
56218672 923 {
924
516043cd 925 const AliAlignObjMatrix* oldAlignObj(0);
926 const Int_t moduleId( kModuleTransformer->GetModuleId() );
927 const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId );
928 for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
56218672 929 {
516043cd 930
931 const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
932 if( localAlignObj && localAlignObj->GetVolUID() == volId )
56218672 933 {
516043cd 934 oldAlignObj = localAlignObj;
935 break;
56218672 936 }
937
010eb601 938 }
56218672 939
516043cd 940 // multiply
941 if( oldAlignObj )
942 {
943
944 TGeoHMatrix oldMatrix;
945 oldAlignObj->GetMatrix( oldMatrix );
946 deltaModuleTransform.Multiply( &oldMatrix );
947
948 }
949
010eb601 950 }
56218672 951
516043cd 952 // Create module mis alignment matrix
953 newGeometryTransformer ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
56218672 954
516043cd 955 AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
010eb601 956
516043cd 957 TIter next(detElements->CreateIterator());
958 AliMUONGeometryDetElement* detElement;
959 Int_t iDe(-1);
960 while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
961 {
962 ++iDe;
963 // make a new detection element
964 AliMUONGeometryDetElement *newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
965 TString lDetElemName(detElement->GetDEName());
966 lDetElemName.ReplaceAll("DE","");
010eb601 967
516043cd 968 // store detector element id and number
969 const Int_t iDetElemId = lDetElemName.Atoi();
970 if( !DetElemIsValid( iDetElemId ) )
971 {
972 AliInfo( Form( "Skipping invalid detector element %i", iDetElemId ) );
973 continue;
974 }
64e2c144 975
516043cd 976 const Int_t iDetElemNumber( GetDetElemNumber( iDetElemId ) );
64e2c144 977
516043cd 978 for( int i=0; i<fgNParCh; ++i )
979 {
980 lDetElemMisAlignment[i] = 0.0;
981 if( iMt<fgNTrkMod ) { lDetElemMisAlignment[i] = misAlignments[iDetElemNumber*fgNParCh+i]; }
982 }
64e2c144 983
516043cd 984 // get transformation
985 TGeoHMatrix deltaGlobalTransform( DeltaTransform( lDetElemMisAlignment ) );
56218672 986
516043cd 987 // update module
988 TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
989 TGeoHMatrix newGlobalTransform( AliMUONGeometryBuilder::Multiply( deltaGlobalTransform, globalTransform ) );
990 newDetElement->SetGlobalTransformation( newGlobalTransform );
991 newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
56218672 992
516043cd 993 // Get matching old alignment and update current matrix accordingly
994 if( oldMisAlignArray )
995 {
56218672 996
516043cd 997 const AliAlignObjMatrix* oldAlignObj(0);
998 const int detElemId( detElement->GetId() );
999 const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
1000 for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
1001 {
56218672 1002
516043cd 1003 const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
1004 if( localAlignObj && localAlignObj->GetVolUID() == volId )
1005 {
1006 oldAlignObj = localAlignObj;
1007 break;
1008 }
56218672 1009
516043cd 1010 }
56218672 1011
516043cd 1012 // multiply
1013 if( oldAlignObj )
1014 {
56218672 1015
516043cd 1016 TGeoHMatrix oldMatrix;
1017 oldAlignObj->GetMatrix( oldMatrix );
1018 deltaGlobalTransform.Multiply( &oldMatrix );
56218672 1019
516043cd 1020 }
56218672 1021
516043cd 1022 }
010eb601 1023
516043cd 1024 // Create misalignment matrix
1025 newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
010eb601 1026
516043cd 1027 }
56218672 1028
516043cd 1029 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1030 }
56218672 1031
516043cd 1032 return newGeometryTransformer;
56218672 1033
010eb601 1034}
1035
56218672 1036//______________________________________________________________________
516043cd 1037void AliMUONAlignment::SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY )
56218672 1038{
1039
516043cd 1040 /// Set alignment resolution to misalign objects to be stored in CDB
1041 /// if rChId is > 0 set parameters for this chamber only, counting from 1
1042 TMatrixDSym mChCorrMatrix(6);
1043 mChCorrMatrix[0][0]=chResX*chResX;
1044 mChCorrMatrix[1][1]=chResY*chResY;
56218672 1045
516043cd 1046 TMatrixDSym mDECorrMatrix(6);
1047 mDECorrMatrix[0][0]=deResX*deResX;
1048 mDECorrMatrix[1][1]=deResY*deResY;
010eb601 1049
516043cd 1050 AliAlignObjMatrix *alignMat = 0x0;
56218672 1051
516043cd 1052 for( Int_t chId = 0; chId <= 9; ++chId )
56218672 1053 {
1054
516043cd 1055 // skip chamber if selection is valid, and does not match
1056 if( rChId > 0 && chId+1 != rChId ) continue;
56218672 1057
516043cd 1058 TString chName1;
1059 TString chName2;
1060 if (chId<4)
1061 {
56218672 1062
516043cd 1063 chName1 = Form("GM%d",chId);
1064 chName2 = Form("GM%d",chId);
56218672 1065
516043cd 1066 } else {
56218672 1067
516043cd 1068 chName1 = Form("GM%d",4+(chId-4)*2);
1069 chName2 = Form("GM%d",4+(chId-4)*2+1);
56218672 1070
516043cd 1071 }
1072
1073 for( int i=0; i<misAlignArray->GetEntries(); ++i )
1074 {
1075
1076 alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
1077 TString volName(alignMat->GetSymName());
1078 if((volName.Contains(chName1)&&
1079 ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
1080 (volName.Length()==volName.Index(chName1)+chName1.Length())))||
1081 (volName.Contains(chName2)&&
1082 ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
1083 (volName.Length()==volName.Index(chName2)+chName2.Length()))))
1084 {
1085
1086 volName.Remove(0,volName.Last('/')+1);
1087 if (volName.Contains("GM")) alignMat->SetCorrMatrix(mChCorrMatrix);
1088 else if (volName.Contains("DE")) alignMat->SetCorrMatrix(mDECorrMatrix);
1089
1090 }
1091
1092 }
1093
1094 }
56218672 1095
1096}
1097
516043cd 1098
56218672 1099//_____________________________________________________
516043cd 1100void AliMUONAlignment::FillDetElemData( AliMUONVCluster* cluster )
56218672 1101{
1102
516043cd 1103 /// Get information of current detection element
1104 // get detector element number from Alice ID
1105 const Int_t detElemId = cluster->GetDetElemId();
1106 fDetElemNumber = GetDetElemNumber( detElemId );
1107
1108 // get detector element
1109 const AliMUONGeometryDetElement* detElement = fTransform->GetDetElement( detElemId );
010eb601 1110
516043cd 1111 /*
1112 get the global transformation matrix and store its inverse, in order to manually perform
1113 the global to Local transformations needed to calculate the derivatives
1114 */
1115 fGeoCombiTransInverse = detElement->GetGlobalTransformation()->Inverse();
56218672 1116
010eb601 1117}
1118
56218672 1119//______________________________________________________________________
516043cd 1120void AliMUONAlignment::FillRecPointData( AliMUONVCluster* cluster )
56218672 1121{
1122
96ebe67e 1123 /// Get information of current cluster
516043cd 1124 fClustPos[0] = cluster->GetX();
1125 fClustPos[1] = cluster->GetY();
1126 fClustPos[2] = cluster->GetZ();
56218672 1127
010eb601 1128}
1129
56218672 1130//______________________________________________________________________
516043cd 1131void AliMUONAlignment::FillTrackParamData( AliMUONTrackParam* trackParam )
56218672 1132{
1133
96ebe67e 1134 /// Get information of current track at current cluster
516043cd 1135 fTrackPos[0] = trackParam->GetNonBendingCoor();
1136 fTrackPos[1] = trackParam->GetBendingCoor();
1137 fTrackPos[2] = trackParam->GetZ();
1138 fTrackSlope[0] = trackParam->GetNonBendingSlope();
1139 fTrackSlope[1] = trackParam->GetBendingSlope();
56218672 1140
010eb601 1141}
1142
516043cd 1143//______________________________________________________________________
516043cd 1144void AliMUONAlignment::LocalEquationX( void )
1145{
1146 /// local equation along X
1147
1148 // 'inverse' (GlobalToLocal) rotation matrix
1149 const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
1150
1151 // local derivatives
1152 SetLocalDerivative( 0, r[0] );
1153 SetLocalDerivative( 1, r[0]*(fTrackPos[2] - fTrackPos0[2]) );
1154 SetLocalDerivative( 2, r[1] );
1155 SetLocalDerivative( 3, r[1]*(fTrackPos[2] - fTrackPos0[2]) );
1156
1157 // global derivatives
1158 /*
1159 alignment parameters are
1160 0: delta_x
1161 1: delta_y
1162 2: delta_phiz
1163 3: delta_z
1164 */
1165
1166 SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[0] );
1167 SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[1] );
1168
1169 if( fBFieldOn )
1170 {
010eb601 1171
516043cd 1172 // use local position for derivatives vs 'delta_phi_z'
1173 SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*fTrackPos[0] + r[0]*fTrackPos[1] );
56218672 1174
516043cd 1175 // use local slopes for derivatives vs 'delta_z'
1176 SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope[0] + r[1]*fTrackSlope[1] );
56218672 1177
516043cd 1178 } else {
56218672 1179
516043cd 1180 // local copy of extrapolated track positions
1181 const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1182 const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
56218672 1183
516043cd 1184 // use properly extrapolated position for derivatives vs 'delta_phi_z'
1185 SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*trackPosX + r[0]*trackPosY );
56218672 1186
516043cd 1187 // use slopes at origin for derivatives vs 'delta_z'
1188 SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[0]*fTrackSlope0[0] + r[1]*fTrackSlope0[1] );
56218672 1189
516043cd 1190 }
56218672 1191
516043cd 1192 // store local equation
1193 fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0] );
56218672 1194
516043cd 1195}
56218672 1196
516043cd 1197//______________________________________________________________________
1198void AliMUONAlignment::LocalEquationY( void )
1199{
1200 /// local equation along Y
56218672 1201
516043cd 1202 // 'inverse' (GlobalToLocal) rotation matrix
1203 const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
56218672 1204
516043cd 1205 // store local derivatives
1206 SetLocalDerivative( 0, r[3] );
1207 SetLocalDerivative( 1, r[3]*(fTrackPos[2] - fTrackPos0[2] ) );
1208 SetLocalDerivative( 2, r[4] );
1209 SetLocalDerivative( 3, r[4]*(fTrackPos[2] - fTrackPos0[2] ) );
56218672 1210
516043cd 1211 // set global derivatives
1212 SetGlobalDerivative( fDetElemNumber*fgNParCh + 0, -r[3]);
1213 SetGlobalDerivative( fDetElemNumber*fgNParCh + 1, -r[4]);
56218672 1214
516043cd 1215 if( fBFieldOn )
1216 {
56218672 1217
516043cd 1218 // use local position for derivatives vs 'delta_phi'
1219 SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*fTrackPos[0] + r[3]*fTrackPos[1]);
56218672 1220
516043cd 1221 // use local slopes for derivatives vs 'delta_z'
1222 SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope[0]+r[4]*fTrackSlope[1] );
010eb601 1223
516043cd 1224 } else {
56218672 1225
516043cd 1226 // local copy of extrapolated track positions
1227 const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1228 const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
56218672 1229
516043cd 1230 // use properly extrapolated position for derivatives vs 'delta_phi'
1231 SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*trackPosX + r[3]*trackPosY );
56218672 1232
516043cd 1233 // use slopes at origin for derivatives vs 'delta_z'
1234 SetGlobalDerivative( fDetElemNumber*fgNParCh + 3, r[3]*fTrackSlope0[0]+r[4]*fTrackSlope0[1] );
56218672 1235
516043cd 1236 }
56218672 1237
516043cd 1238 // store local equation
1239 fMillepede->SetLocalEquation( fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1] );
56218672 1240
010eb601 1241}
1242
516043cd 1243//_________________________________________________________________________
1244TGeoCombiTrans AliMUONAlignment::DeltaTransform( const double *lMisAlignment) const
56218672 1245{
516043cd 1246 /// Get Delta Transformation, based on alignment parameters
56218672 1247
516043cd 1248 // translation
1249 const TGeoTranslation deltaTrans( lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
3b2890be 1250
516043cd 1251 // rotation
1252 TGeoRotation deltaRot;
1253 deltaRot.RotateZ(lMisAlignment[2]*180./TMath::Pi());
56218672 1254
516043cd 1255 // combined rotation and translation.
1256 return TGeoCombiTrans(deltaTrans,deltaRot);
56218672 1257
010eb601 1258}
1259
516043cd 1260//______________________________________________________________________
1261void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value)
56218672 1262{
516043cd 1263 /// Constrain equation defined by par to value
1264 if( !fInitialized )
1265 { AliFatal( "Millepede is not initialized" ); }
010eb601 1266
516043cd 1267 fMillepede->SetGlobalConstraint(par, value);
010eb601 1268}
1269
516043cd 1270//______________________________________________________________________
1271Bool_t AliMUONAlignment::DetElemIsValid( Int_t iDetElemId ) const
010eb601 1272{
516043cd 1273 /// return true if given detector element is valid (and belongs to muon tracker)
1274 const Int_t iCh = iDetElemId/100;
1275 const Int_t iDet = iDetElemId%100;
1276 return ( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] );
010eb601 1277}
1278
1279//______________________________________________________________________
516043cd 1280Int_t AliMUONAlignment::GetDetElemNumber( Int_t iDetElemId ) const
010eb601 1281{
516043cd 1282 /// get det element number from ID
1283 // get chamber and element number in chamber
1284 const Int_t iCh = iDetElemId/100;
1285 const Int_t iDet = iDetElemId%100;
56218672 1286
516043cd 1287 // make sure detector index is valid
1288 if( !( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] ) )
1289 { AliFatal( Form( "Invalid detector element id: %i", iDetElemId ) ); }
3b2890be 1290
516043cd 1291 // add number of detectors up to this chamber
1292 return iDet + fgSNDetElemCh[iCh-1];
56218672 1293
010eb601 1294}
1295
4818a9b7 1296//______________________________________________________________________
516043cd 1297Int_t AliMUONAlignment::GetChamberId( Int_t iDetElemNumber ) const
56218672 1298{
516043cd 1299 /// get chamber (counting from 1) matching a given detector element id
1300 Int_t iCh( 0 );
1301 for( iCh=0; iCh<fgNCh; iCh++ )
1302 { if( iDetElemNumber < fgSNDetElemCh[iCh] ) break; }
56218672 1303
516043cd 1304 return iCh;
1305}
4818a9b7 1306
516043cd 1307//______________________________________________________________________
1308TString AliMUONAlignment::GetParameterMaskString( UInt_t mask ) const
1309{
1310 TString out;
1311 if( mask & ParX ) out += "X";
1312 if( mask & ParY ) out += "Y";
1313 if( mask & ParZ ) out += "Z";
1314 if( mask & ParTZ ) out += "T";
1315 return out;
1316}
56218672 1317
516043cd 1318//______________________________________________________________________
1319TString AliMUONAlignment::GetSidesMaskString( UInt_t mask ) const
1320{
1321 TString out;
1322 if( mask & SideTop ) out += "T";
1323 if( mask & SideLeft ) out += "L";
1324 if( mask & SideBottom ) out += "B";
1325 if( mask & SideRight ) out += "R";
1326 return out;
4818a9b7 1327}