AliRoot Core  edcc906 (edcc906)
AliMUONAlignment.cxx
Go to the documentation of this file.
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: AliMUONAlignment.cxx 51000 2011-08-08 17:58:17Z ivana $ */
17 
18 //-----------------------------------------------------------------------------
28 //-----------------------------------------------------------------------------
29 
30 #include "AliMUONAlignment.h"
31 #include "AliMUONTrack.h"
32 #include "AliMUONTrackParam.h"
33 #include "AliMUONVCluster.h"
37 #include "AliMUONGeometryBuilder.h"
38 #include "AliMillePede2.h"
39 
40 #include "AliMpExMap.h"
41 #include "AliMpExMapIterator.h"
42 
43 #include "AliAlignObjMatrix.h"
44 #include "AliLog.h"
45 
46 #include <TMath.h>
47 #include <TMatrixDSym.h>
48 #include <TMatrixD.h>
49 #include <TClonesArray.h>
50 #include <TGraphErrors.h>
51 
53 ClassImp(AliMUONAlignment)
54 ClassImp(LocalTrackParam)
56 
57 //_____________________________________________________________________
58 // static variables
59 const Int_t AliMUONAlignment::fgNDetElemCh[AliMUONAlignment::fgNCh] = { 4, 4, 4, 4, 18, 18, 26, 26, 26, 26 };
60 const Int_t AliMUONAlignment::fgSNDetElemCh[AliMUONAlignment::fgNCh+1] = { 0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156 };
61 
62 
63 // number of detector elements in each half-chamber
64 const Int_t AliMUONAlignment::fgNDetElemHalfCh[AliMUONAlignment::fgNHalfCh] = { 2, 2, 2, 2, 2, 2, 2, 2, 9, 9, 9, 9, 13, 13, 13, 13, 13, 13, 13, 13 };
65 
66 // list of detector elements for each half chamber
68 {
69 { 100, 103, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
70 { 101, 102, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
71 
72 { 200, 203, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
73 { 201, 202, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
74 
75 { 300, 303, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
76 { 301, 302, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
77 
78 { 400, 403, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
79 { 401, 402, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
80 
81 { 500, 501, 502, 503, 504, 514, 515, 516, 517, 0, 0, 0, 0 },
82 { 505, 506, 507, 508, 509, 510, 511, 512, 513, 0, 0, 0, 0 },
83 
84 { 600, 601, 602, 603, 604, 614, 615, 616, 617, 0, 0, 0, 0 },
85 { 605, 606, 607, 608, 609, 610, 611, 612, 613, 0, 0, 0, 0 },
86 
87 { 700, 701, 702, 703, 704, 705, 706, 720, 721, 722, 723, 724, 725 },
88 { 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719 },
89 
90 { 800, 801, 802, 803, 804, 805, 806, 820, 821, 822, 823, 824, 825 },
91 { 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819 },
92 
93 { 900, 901, 902, 903, 904, 905, 906, 920, 921, 922, 923, 924, 925 },
94 { 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919 },
95 
96 { 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1020, 1021, 1022, 1023, 1024, 1025 },
97 { 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019 }
98 
99 };
100 
101 //_____________________________________________________________________
103 class Array
104 {
105 
106  public:
107 
109  Array( void )
110  {
111  for( Int_t i=0; i < AliMUONAlignment::fNGlobal; ++i )
112  { values[i] = 0; }
113  }
114 
117 
118  private:
119 
121  Array(const Array& );
122 
124  Array& operator = (const Array& );
125 
126 };
127 
128 //________________________________________________________________________
129 Double_t Square( Double_t x ) { return x*x; }
130 
131 //_____________________________________________________________________
133  : TObject(),
134  fInitialized( kFALSE ),
135  fRunNumber( 0 ),
136  fBFieldOn( kFALSE ),
137  fRefitStraightTracks( kFALSE ),
138  fStartFac( 256 ),
139  fResCutInitial( 100 ),
140  fResCut( 100 ),
141  fMillepede( 0L ),
142  fCluster( 0L ),
143  fNStdDev( 3 ),
144  fDetElemNumber( 0 ),
145  fTrackRecord(),
146  fTransform( 0 ),
147  fGeoCombiTransInverse(),
148  fDoEvaluation( kFALSE ),
149  fTrackParamOrig( 0 ),
150  fTrackParamNew( 0 ),
151  fTFile( 0 ),
152  fTTree( 0 )
153 {
155  fSigma[0] = 1.5e-1;
156  fSigma[1] = 1.0e-2;
157 
158  // default allowed variations
159  fAllowVar[0] = 0.5; // x
160  fAllowVar[1] = 0.5; // y
161  fAllowVar[2] = 0.01; // phi_z
162  fAllowVar[3] = 5; // z
163 
164  // initialize millepede
165  fMillepede = new AliMillePede2();
166 
167  // initialize degrees of freedom
168  // by default all parameters are free
169  for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
171 
172  // initialize local equations
173  for(int i=0; i<fNLocal; ++i )
174  { fLocalDerivatives[i] = 0.0; }
175 
176  for(int i=0; i<fNGlobal; ++i )
177  { fGlobalDerivatives[i] = 0.0; }
178 
179 }
180 
181 //_____________________________________________________________________
183 {
185 }
186 
187 //_____________________________________________________________________
189 {
190 
192 
197  if( fInitialized )
198  { AliFatal( "Millepede already initialized" ); }
199 
200  // assign proper groupID to free parameters
201  Int_t nGlobal = 0;
202  for( Int_t iPar = 0; iPar < fNGlobal; ++iPar )
203  {
204 
205  if( fGlobalParameterStatus[iPar] == kFixedParId )
206  {
207  // fixed parameters are left unchanged
208  continue;
209 
210  } else if( fGlobalParameterStatus[iPar] == kFreeParId || fGlobalParameterStatus[iPar] == kGroupBaseId ) {
211 
212  // free parameters or first element of group are assigned a new group id
213  fGlobalParameterStatus[iPar] = nGlobal++;
214  continue;
215 
216  } else if( fGlobalParameterStatus[iPar] < kGroupBaseId ) {
217 
218  // get detector element id from status, get chamber parameter id
219  const Int_t iDeBase( kGroupBaseId - 1 - fGlobalParameterStatus[iPar] );
220  const Int_t iParBase = iPar%fgNParCh;
221 
222  // check
223  if( iDeBase < 0 || iDeBase >= iPar/fgNParCh )
224  { AliFatal( Form( "Group for parameter index %i has wrong base detector element: %i", iPar, iDeBase ) ); }
225 
226  // assign identical group id to current
227  fGlobalParameterStatus[iPar] = fGlobalParameterStatus[iDeBase*fgNParCh + iParBase];
228  AliInfo( Form( "Parameter %i grouped to detector %i (%s)", iPar, iDeBase, GetParameterMaskString( 1<<iParBase ).Data() ) );
229 
230  } else AliFatal( Form( "Unrecognized parameter status for index %i: %i", iPar, fGlobalParameterStatus[iPar] ) );
231 
232  }
233 
234  AliInfo( Form( "Free Parameters: %i out of %i", nGlobal, fNGlobal ) );
235 
236  // initialize millepede
238  fInitialized = kTRUE;
239 
240  // some debug output
241  for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
242  { AliInfo( Form( "fAllowVar[%i]= %f", iPar, fAllowVar[iPar] ) ); }
243 
244  // set allowed variations for all parameters
245  for( Int_t iDet = 0; iDet < fgNDetElem; ++iDet )
246  {
247  for( Int_t iPar = 0; iPar < fgNParCh; ++iPar )
248  { fMillepede->SetParSigma( iDet*fgNParCh + iPar, fAllowVar[iPar] ); }
249  }
250 
251  // Set iterations
253 
254  // setup monitoring TFile
256  {
257  fTFile = new TFile( "AliMUONAlignment.root","RECREATE");
258  fTTree = new TTree("TreeE","Evaluation");
259 
260  const Int_t kSplitlevel = 98;
261  const Int_t kBufsize = 32000;
262 
264  fTTree->Branch( "fTrackParamOrig", "LocalTrackParam", &fTrackParamOrig, kBufsize, kSplitlevel );
265 
267  fTTree->Branch( "fTrackParamNew", "LocalTrackParam", &fTrackParamNew, kBufsize, kSplitlevel );
268  }
269 
270 }
271 
272 //_____________________________________________________
274 {
275  AliInfo( "Closing Evaluation TFile" );
276  if( fTFile && fTTree )
277  {
278  fTFile->cd();
279  fTTree->Write();
280  fTFile->Close();
281  }
282 }
283 
284 //_____________________________________________________
285 AliMillePedeRecord* AliMUONAlignment::ProcessTrack( AliMUONTrack* track, Bool_t doAlignment, Double_t weight )
286 {
288 
293  // reset track records
296 
297  // loop over clusters to get starting values
298  Bool_t first( kTRUE );
299  TObjArray* trackParamAtCluster( track->GetTrackParamAtCluster() );
300  for( Int_t iTrackParam = 0; iTrackParam < trackParamAtCluster->GetEntries(); ++iTrackParam )
301  {
302 
303  // get track parameters
304  AliMUONTrackParam* trackParam( (AliMUONTrackParam *) trackParamAtCluster->At(iTrackParam) );
305  if( !trackParam ) continue;
306 
307  // get cluster
308  AliMUONVCluster* cluster = trackParam->GetClusterPtr();
309  if( !cluster ) continue;
310 
311  // for first valid cluster, save track position as "starting" values
312  if( first )
313  {
314 
315  first = kFALSE;
316  FillTrackParamData( trackParam );
317  fTrackPos0[0] = fTrackPos[0];
318  fTrackPos0[1] = fTrackPos[1];
319  fTrackPos0[2] = fTrackPos[2];
320  fTrackSlope0[0] = fTrackSlope[0];
321  fTrackSlope0[1] = fTrackSlope[1];
322 
323  break;
324 
325  }
326 
327  }
328 
329  // redo straight track fit
331  {
332 
333  // refit straight track
334  const LocalTrackParam trackParam( RefitStraightTrack( track, fTrackPos0[2] ) );
335 
336  // fill evaluation tree
337  if( fTrackParamOrig )
338  {
344  }
345 
346  // new ones
347  if( fTrackParamNew )
348  {
349  fTrackParamNew->fTrackX = trackParam.fTrackX;
350  fTrackParamNew->fTrackY = trackParam.fTrackY;
351  fTrackParamNew->fTrackZ = trackParam.fTrackZ;
354  }
355 
356  if( fTTree ) fTTree->Fill();
357 
358  /*
359  copy new parameters to stored ones for derivatives calculation
360  this is done only if BFieldOn is false, for which these parameters are used
361  */
362  if( !fBFieldOn )
363  {
364  fTrackPos0[0] = trackParam.fTrackX;
365  fTrackPos0[1] = trackParam.fTrackY;
366  fTrackPos0[2] = trackParam.fTrackZ;
367  fTrackSlope[0] = trackParam.fTrackSlopeX;
368  fTrackSlope[1] = trackParam.fTrackSlopeY;
369  }
370  }
371 
372  // second loop to perform alignment
373  for( Int_t iTrackParam = 0; iTrackParam < trackParamAtCluster->GetEntries(); ++iTrackParam )
374  {
375 
376  // get track parameters
377  AliMUONTrackParam* trackParam( (AliMUONTrackParam *) trackParamAtCluster->At(iTrackParam) );
378  if( !trackParam ) continue;
379 
380  // get cluster
381  AliMUONVCluster* cluster = trackParam->GetClusterPtr();
382  if( !cluster ) continue;
383 
384  // fill local variables for this position --> one measurement
385  FillDetElemData( cluster );
386  FillRecPointData( cluster );
387  FillTrackParamData( trackParam );
388 
389  // 'inverse' (GlobalToLocal) rotation matrix
390  const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
391 
392  // calculate measurements
393  if( fBFieldOn )
394  {
395 
396  // use residuals (cluster - track) for measurement
397  fMeas[0] = r[0]*(fClustPos[0] - fTrackPos[0]) + r[1]*(fClustPos[1] - fTrackPos[1]);
398  fMeas[1] = r[3]*(fClustPos[0] - fTrackPos[0]) + r[4]*(fClustPos[1] - fTrackPos[1]);
399 
400  } else {
401 
402  // use cluster position for measurement
403  fMeas[0] = ( r[0]*fClustPos[0] + r[1]*fClustPos[1] );
404  fMeas[1] = ( r[3]*fClustPos[0] + r[4]*fClustPos[1] );
405 
406  }
407 
408  // Set local equations
409  LocalEquationX();
410  LocalEquationY();
411 
412  }
413 
414  // copy track record
416  fMillepede->SetRecordWeight(weight);
418 
419  // save record data
420  if( doAlignment )
421  {
423  }
424 
425  // return record
426  return &fTrackRecord;
427 
428 }
429 
430 //______________________________________________________________________________
432 {
434  if( !trackRecord ) return;
435 
436  // make sure record storage is initialized
438 
439  // copy content
440  *fMillepede->GetRecord() = *trackRecord;
441 
442  // save record
444 
445  return;
446 
447 }
448 
449 //_____________________________________________________________________
450 void AliMUONAlignment::FixAll( UInt_t mask )
451 {
453  AliInfo( Form( "Fixing %s for all detector elements", GetParameterMaskString( mask ).Data() ) );
454 
455  // fix all stations
456  for( Int_t i = 0; i < fgNDetElem; ++i )
457  {
458  if( mask & ParX ) FixParameter(i, 0);
459  if( mask & ParY ) FixParameter(i, 1);
460  if( mask & ParTZ ) FixParameter(i, 2);
461  if( mask & ParZ ) FixParameter(i, 3);
462  }
463 
464 }
465 
466 //_____________________________________________________________________
467 void AliMUONAlignment::FixChamber( Int_t iCh, UInt_t mask )
468 {
470 
471  // check boundaries
472  if( iCh < 1 || iCh > 10 )
473  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
474 
475  // get first and last element
476  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
477  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
478  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
479  {
480 
481  AliInfo( Form( "Fixing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
482 
483  if( mask & ParX ) FixParameter(i, 0);
484  if( mask & ParY ) FixParameter(i, 1);
485  if( mask & ParTZ ) FixParameter(i, 2);
486  if( mask & ParZ ) FixParameter(i, 3);
487 
488  }
489 }
490 
491 //_____________________________________________________________________
492 void AliMUONAlignment::FixDetElem( Int_t iDetElemId, UInt_t mask )
493 {
495  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
496  if ( mask & ParX ) FixParameter(iDet, 0);
497  if ( mask & ParY ) FixParameter(iDet, 1);
498  if ( mask & ParTZ ) FixParameter(iDet, 2);
499  if ( mask & ParZ ) FixParameter(iDet, 3);
500 
501 }
502 
503 //_____________________________________________________________________
504 void AliMUONAlignment::FixHalfSpectrometer( const Bool_t *lChOnOff, UInt_t sidesMask, UInt_t mask )
505 {
506 
508  for( Int_t i = 0; i < fgNDetElem; ++i )
509  {
510 
511  // get chamber matching detector
512  const Int_t iCh( GetChamberId(i) );
513  if( !lChOnOff[iCh-1] ) continue;
514 
515  // get detector element in chamber
516  Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
517 
518  // skip detector if its side is off
519  // stations 1 and 2
520  if( iCh>=1 && iCh<=4 )
521  {
522  if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
523  if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
524  if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
525  if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
526  }
527 
528  // station 3
529  if (iCh>=5 && iCh<=6)
530  {
531  if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
532  if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
533  if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
534  if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
535  }
536 
537  // stations 4 and 5
538  if (iCh>=7 && iCh<=10)
539  {
540  if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
541  if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
542  if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
543  if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
544  }
545 
546  // detector is accepted, fix it
547  FixDetElem( i, mask );
548 
549  }
550 
551 }
552 
553 //______________________________________________________________________
555 {
556 
558  if( fInitialized )
559  { AliFatal( "Millepede already initialized" ); }
560 
562 
563 }
564 
565 
566 //_____________________________________________________________________
567 void AliMUONAlignment::ReleaseChamber( Int_t iCh, UInt_t mask )
568 {
570 
571  // check boundaries
572  if( iCh < 1 || iCh > 10 )
573  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
574 
575  // get first and last element
576  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
577  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
578  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
579  {
580 
581  AliInfo( Form( "Releasing %s for detector element %i", GetParameterMaskString(mask).Data(), i ) );
582 
583  if( mask & ParX ) ReleaseParameter(i, 0);
584  if( mask & ParY ) ReleaseParameter(i, 1);
585  if( mask & ParTZ ) ReleaseParameter(i, 2);
586  if( mask & ParZ ) ReleaseParameter(i, 3);
587 
588  }
589 }
590 
591 //_____________________________________________________________________
592 void AliMUONAlignment::ReleaseDetElem( Int_t iDetElemId, UInt_t mask )
593 {
595  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
596  if ( mask & ParX ) ReleaseParameter(iDet, 0);
597  if ( mask & ParY ) ReleaseParameter(iDet, 1);
598  if ( mask & ParTZ ) ReleaseParameter(iDet, 2);
599  if ( mask & ParZ ) ReleaseParameter(iDet, 3);
600 
601 }
602 
603 //______________________________________________________________________
605 {
606 
608  if( fInitialized )
609  { AliFatal( "Millepede already initialized" ); }
610 
612 
613 }
614 
615 //_____________________________________________________________________
616 void AliMUONAlignment::GroupChamber( Int_t iCh, UInt_t mask )
617 {
619  if( iCh < 1 || iCh > fgNCh )
620  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
621 
622  const Int_t detElemMin = 100*iCh;
623  const Int_t detElemMax = 100*iCh + fgNDetElemCh[iCh]-1;
624  GroupDetElems( detElemMin, detElemMax, mask );
625 
626 }
627 
628 //_____________________________________________________________________
629 void AliMUONAlignment::GroupHalfChamber( Int_t iCh, Int_t iHalf, UInt_t mask )
630 {
632  if( iCh < 1 || iCh > fgNCh )
633  { AliFatal( Form( "Invalid chamber index %i", iCh ) ); }
634 
635  if( iHalf < 0 || iHalf > 1 )
636  { AliFatal( Form( "Invalid half chamber index %i", iHalf ) ); }
637 
638  const Int_t iHalfCh = 2*(iCh-1)+iHalf;
639  GroupDetElems( &fgDetElemHalfCh[iHalfCh][0], fgNDetElemHalfCh[iHalfCh], mask );
640 
641 }
642 
643 //_____________________________________________________________________
644 void AliMUONAlignment::GroupDetElems( Int_t detElemMin, Int_t detElemMax, UInt_t mask )
645 {
647  // check number of detector elements
648  const Int_t nDetElem = detElemMax - detElemMin + 1;
649  if( nDetElem<2 )
650  { AliFatal( Form( "Requested group of DEs %d-%d contains less than 2 DE's", detElemMin, detElemMax ) ); }
651 
652  // create list
653  Int_t* detElemList = new int[nDetElem];
654  for( Int_t i = 0; i < nDetElem; ++i )
655  { detElemList[i] = detElemMin+i; }
656 
657  // group
658  GroupDetElems( detElemList, nDetElem, mask );
659  delete[] detElemList;
660 
661 }
662 
663 //_____________________________________________________________________
664 void AliMUONAlignment::GroupDetElems( const Int_t* detElemList, Int_t nDetElem, UInt_t mask )
665 {
667  if( fInitialized )
668  { AliFatal( "Millepede already initialized" ); }
669 
670  const Int_t iDeBase( GetDetElemNumber( detElemList[0] ) );
671  for( Int_t i = 0; i < nDetElem; ++i )
672  {
673  const Int_t iDeCurrent( GetDetElemNumber( detElemList[i] ) );
674  if( mask & ParX ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 0] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
675  if( mask & ParY ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 1] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
676  if( mask & ParTZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 2] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
677  if( mask & ParZ ) fGlobalParameterStatus[iDeCurrent*fgNParCh + 3] = (i==0) ? kGroupBaseId : (kGroupBaseId-iDeBase-1);
678 
679  if( i== 0 ) AliInfo( Form( "Creating new group for detector %i and variable %s", detElemList[i], GetParameterMaskString( mask ).Data() ) );
680  else AliInfo( Form( "Adding detector element %i to current group", detElemList[i] ) );
681  }
682 
683 }
684 
685 //______________________________________________________________________
686 void AliMUONAlignment::SetChamberNonLinear( Int_t iCh, UInt_t mask )
687 {
689  const Int_t iDetElemFirst = fgSNDetElemCh[iCh-1];
690  const Int_t iDetElemLast = fgSNDetElemCh[iCh];
691  for( Int_t i = iDetElemFirst; i < iDetElemLast; ++i )
692  {
693 
694  if( mask & ParX ) SetParameterNonLinear(i, 0);
695  if( mask & ParY ) SetParameterNonLinear(i, 1);
696  if( mask & ParTZ ) SetParameterNonLinear(i, 2);
697  if( mask & ParZ ) SetParameterNonLinear(i, 3);
698 
699  }
700 
701 }
702 
703 //_____________________________________________________________________
704 void AliMUONAlignment::SetDetElemNonLinear( Int_t iDetElemId, UInt_t mask )
705 {
707  const Int_t iDet( GetDetElemNumber( iDetElemId ) );
708  if ( mask & ParX ) SetParameterNonLinear(iDet, 0);
709  if ( mask & ParY ) SetParameterNonLinear(iDet, 1);
710  if ( mask & ParTZ ) SetParameterNonLinear(iDet, 2);
711  if ( mask & ParZ ) SetParameterNonLinear(iDet, 3);
712 
713 }
714 
715 //______________________________________________________________________
717 {
719  if( !fInitialized )
720  { AliFatal( "Millepede not initialized" ); }
721 
722  fMillepede->SetNonLinear( iPar );
723  AliInfo( Form( "Parameter %i set to non linear", iPar ) );
724 }
725 
726 //______________________________________________________________________
727 void AliMUONAlignment::AddConstraints( const Bool_t *lChOnOff, UInt_t mask )
728 {
730 
731  Array fConstraintX;
732  Array fConstraintY;
733  Array fConstraintTZ;
734  Array fConstraintZ;
735 
736  for( Int_t i = 0; i < fgNDetElem; ++i )
737  {
738 
739  // get chamber matching detector
740  const Int_t iCh( GetChamberId(i) );
741  if (lChOnOff[iCh-1])
742  {
743 
744  if( mask & ParX ) fConstraintX.values[i*fgNParCh+0]=1.0;
745  if( mask & ParY ) fConstraintY.values[i*fgNParCh+1]=1.0;
746  if( mask & ParTZ ) fConstraintTZ.values[i*fgNParCh+2]=1.0;
747  if( mask & ParZ ) fConstraintTZ.values[i*fgNParCh+3]=1.0;
748 
749  }
750  }
751 
752  if( mask & ParX ) AddConstraint(fConstraintX.values,0.0);
753  if( mask & ParY ) AddConstraint(fConstraintY.values,0.0);
754  if( mask & ParTZ ) AddConstraint(fConstraintTZ.values,0.0);
755  if( mask & ParZ ) AddConstraint(fConstraintZ.values,0.0);
756 
757 }
758 
759 //______________________________________________________________________
760 void AliMUONAlignment::AddConstraints(const Bool_t *lChOnOff, const Bool_t *lVarXYT, UInt_t sidesMask )
761 {
762  /*
763  questions:
764  - is there not redundancy/inconsistency between lDetTLBR and lSpecLROnOff ? shouldn't we use only lDetTLBR ?
765  - why is weight ignored for ConstrainT and ConstrainB
766  - why is there no constrain on z
767  */
768 
770  Double_t lMeanY = 0.;
771  Double_t lSigmaY = 0.;
772  Double_t lMeanZ = 0.;
773  Double_t lSigmaZ = 0.;
774  Int_t lNDetElem = 0;
775 
776  for( Int_t i = 0; i < fgNDetElem; ++i )
777  {
778 
779  // get chamber matching detector
780  const Int_t iCh( GetChamberId(i) );
781 
782  // skip detector if chamber is off
783  if( lChOnOff[iCh-1] ) continue;
784 
785  // get detector element id from detector element number
786  const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
787  const Int_t lDetElemId = iCh*100+lDetElemNumber;
788 
789  // skip detector if its side is off
790  // stations 1 and 2
791  if( iCh>=1 && iCh<=4 )
792  {
793  if( lDetElemNumber == 0 && !( sidesMask & SideTopRight ) ) continue;
794  if( lDetElemNumber == 1 && !( sidesMask & SideTopLeft ) ) continue;
795  if( lDetElemNumber == 2 && !( sidesMask & SideBottomLeft ) ) continue;
796  if( lDetElemNumber == 3 && !( sidesMask & SideBottomRight ) ) continue;
797  }
798 
799  // station 3
800  if (iCh>=5 && iCh<=6)
801  {
802  if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !( sidesMask & SideTopRight ) ) continue;
803  if( lDetElemNumber >= 5 && lDetElemNumber <= 10 && !( sidesMask & SideTopLeft ) ) continue;
804  if( lDetElemNumber >= 11 && lDetElemNumber <= 13 && !( sidesMask & SideBottomLeft ) ) continue;
805  if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !( sidesMask & SideBottomRight ) ) continue;
806  }
807 
808  // stations 4 and 5
809  if (iCh>=7 && iCh<=10)
810  {
811  if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !( sidesMask & SideTopRight ) ) continue;
812  if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !( sidesMask & SideTopLeft ) ) continue;
813  if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !( sidesMask & SideBottomLeft ) ) continue;
814  if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !( sidesMask & SideBottomRight ) ) continue;
815  }
816 
817  // get global x, y and z position
818  Double_t lDetElemGloX = 0.;
819  Double_t lDetElemGloY = 0.;
820  Double_t lDetElemGloZ = 0.;
821  fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
822 
823  // increment mean Y, mean Z, sigmas and number of accepted detectors
824  lMeanY += lDetElemGloY;
825  lSigmaY += lDetElemGloY*lDetElemGloY;
826  lMeanZ += lDetElemGloZ;
827  lSigmaZ += lDetElemGloZ*lDetElemGloZ;
828  lNDetElem++;
829 
830  }
831 
832  // calculate mean values
833  lMeanY /= lNDetElem;
834  lSigmaY /= lNDetElem;
835  lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
836  lMeanZ /= lNDetElem;
837  lSigmaZ /= lNDetElem;
838  lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
839  AliInfo( Form( "Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ ) );
840 
841  // create all possible arrays
842  Array fConstraintX[4]; //Array for constraint equation X
843  Array fConstraintY[4]; //Array for constraint equation Y
844  Array fConstraintP[4]; //Array for constraint equation P
845  Array fConstraintXZ[4]; //Array for constraint equation X vs Z
846  Array fConstraintYZ[4]; //Array for constraint equation Y vs Z
847  Array fConstraintPZ[4]; //Array for constraint equation P vs Z
848 
849  // do we really need these ?
850  Array fConstraintXY[4]; //Array for constraint equation X vs Y
851  Array fConstraintYY[4]; //Array for constraint equation Y vs Y
852  Array fConstraintPY[4]; //Array for constraint equation P vs Y
853 
854  // fill Bool_t sides array based on masks, for convenience
855  Bool_t lDetTLBR[4];
856  lDetTLBR[0] = sidesMask & SideTop;
857  lDetTLBR[1] = sidesMask & SideLeft;
858  lDetTLBR[2] = sidesMask & SideBottom;
859  lDetTLBR[3] = sidesMask & SideRight;
860 
861  for( Int_t i = 0; i < fgNDetElem; ++i )
862  {
863 
864  // get chamber matching detector
865  const Int_t iCh( GetChamberId(i) );
866 
867  // skip detector if chamber is off
868  if( !lChOnOff[iCh-1] ) continue;
869 
870  // get detector element id from detector element number
871  const Int_t lDetElemNumber = i-fgSNDetElemCh[iCh-1];
872  const Int_t lDetElemId = iCh*100+lDetElemNumber;
873 
874  // get global x, y and z position
875  Double_t lDetElemGloX = 0.;
876  Double_t lDetElemGloY = 0.;
877  Double_t lDetElemGloZ = 0.;
878  fTransform->Local2Global( lDetElemId, 0, 0, 0, lDetElemGloX, lDetElemGloY, lDetElemGloZ );
879 
880  // loop over sides
881  for( Int_t iSide = 0; iSide < 4; iSide++ )
882  {
883 
884  // skip if side is not selected
885  if( !lDetTLBR[iSide] ) continue;
886 
887  // skip detector if it is not in the selected side
888  // stations 1 and 2
889  if( iCh>=1 && iCh<=4 )
890  {
891  if( lDetElemNumber == 0 && !(iSide == 0 || iSide == 3) ) continue; // top-right
892  if( lDetElemNumber == 1 && !(iSide == 0 || iSide == 1) ) continue; // top-left
893  if( lDetElemNumber == 2 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
894  if( lDetElemNumber == 3 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
895  }
896 
897  // station 3
898  if (iCh>=5 && iCh<=6)
899  {
900  if( lDetElemNumber >= 0 && lDetElemNumber <= 4 && !(iSide == 0 || iSide == 3) ) continue; // top-right
901  if( lDetElemNumber >= 5 && lDetElemNumber <= 9 && !(iSide == 0 || iSide == 1) ) continue; // top-left
902  if( lDetElemNumber >= 10 && lDetElemNumber <= 13 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
903  if( lDetElemNumber >= 14 && lDetElemNumber <= 17 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
904  }
905 
906  // stations 4 and 5
907  if (iCh>=7 && iCh<=10)
908  {
909  if( lDetElemNumber >= 0 && lDetElemNumber <= 6 && !(iSide == 0 || iSide == 3) ) continue; // top-right
910  if( lDetElemNumber >= 7 && lDetElemNumber <= 13 && !(iSide == 0 || iSide == 1) ) continue; // top-left
911  if( lDetElemNumber >= 14 && lDetElemNumber <= 19 && !(iSide == 2 || iSide == 1) ) continue; // bottom-left
912  if( lDetElemNumber >= 20 && lDetElemNumber <= 25 && !(iSide == 2 || iSide == 3) ) continue; // bottom-right
913  }
914 
915  // constrain x
916  if( lVarXYT[0] ) fConstraintX[iSide].values[i*fgNParCh+0] = 1;
917 
918  // constrain y
919  if( lVarXYT[1] ) fConstraintY[iSide].values[i*fgNParCh+1] = 1;
920 
921  // constrain phi (rotation around z)
922  if( lVarXYT[2] ) fConstraintP[iSide].values[i*fgNParCh+2] = 1;
923 
924  // x-z shearing
925  if( lVarXYT[3] ) fConstraintXZ[iSide].values[i*fgNParCh+0] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
926 
927  // y-z shearing
928  if( lVarXYT[4] ) fConstraintYZ[iSide].values[i*fgNParCh+1] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
929 
930  // phi-z shearing
931  if( lVarXYT[5] ) fConstraintPZ[iSide].values[i*fgNParCh+2] = (lDetElemGloZ-lMeanZ)/lSigmaZ;
932 
933  // x-y shearing
934  if( lVarXYT[6] ) fConstraintXY[iSide].values[i*fgNParCh+0] = (lDetElemGloY-lMeanY)/lSigmaY;
935 
936  // y-y shearing
937  if( lVarXYT[7] ) fConstraintYY[iSide].values[i*fgNParCh+1] = (lDetElemGloY-lMeanY)/lSigmaY;
938 
939  // phi-y shearing
940  if( lVarXYT[8] ) fConstraintPY[iSide].values[i*fgNParCh+2] = (lDetElemGloY-lMeanY)/lSigmaY;
941 
942  }
943 
944  }
945 
946  // pass constraints to millepede
947  for( Int_t iSide = 0; iSide < 4; iSide++ )
948  {
949  // skip if side is not selected
950  if( !lDetTLBR[iSide] ) continue;
951 
952  if( lVarXYT[0] ) AddConstraint(fConstraintX[iSide].values,0.0);
953  if( lVarXYT[1] ) AddConstraint(fConstraintY[iSide].values,0.0);
954  if( lVarXYT[2] ) AddConstraint(fConstraintP[iSide].values,0.0);
955  if( lVarXYT[3] ) AddConstraint(fConstraintXZ[iSide].values,0.0);
956  if( lVarXYT[4] ) AddConstraint(fConstraintYZ[iSide].values,0.0);
957  if( lVarXYT[5] ) AddConstraint(fConstraintPZ[iSide].values,0.0);
958  if( lVarXYT[6] ) AddConstraint(fConstraintXY[iSide].values,0.0);
959  if( lVarXYT[7] ) AddConstraint(fConstraintYY[iSide].values,0.0);
960  if( lVarXYT[8] ) AddConstraint(fConstraintPY[iSide].values,0.0);
961  }
962 
963 }
964 
965 //______________________________________________________________________
967 {
969  if( !fInitialized )
970  { AliFatal( "Millepede is not initialized" ); }
971 
973 }
974 
975 //______________________________________________________________________
976 void AliMUONAlignment::SetAllowedVariation( Int_t iPar, Double_t value )
977 {
979  // check initialization
980  if( fInitialized )
981  { AliFatal( "Millepede already initialized" ); }
982 
983  // check initialization
984  if( !(iPar >= 0 && iPar < fgNParCh ) )
985  { AliFatal( Form( "Invalid index: %i", iPar ) ); }
986 
987  fAllowVar[iPar] = value;
988 }
989 
990 //______________________________________________________________________
991 void AliMUONAlignment::SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
992 {
993 
995  fSigma[0] = sigmaX;
996  fSigma[1] = sigmaY;
997 
998  // print
999  for( Int_t i=0; i<2; ++i )
1000  { AliInfo( Form( "fSigma[%i]=%f", i, fSigma[i] ) ); }
1001 
1002 }
1003 
1004 //_____________________________________________________
1005 void AliMUONAlignment::GlobalFit( Double_t *parameters, Double_t *errors, Double_t *pulls )
1006 {
1007 
1009  fMillepede->GlobalFit( parameters, errors, pulls );
1010 
1011  AliInfo( "Done fitting global parameters" );
1012  for( int iDet=0; iDet<fgNDetElem; ++iDet )
1013  {
1014  AliInfo( Form( "%d\t %f\t %f\t %f\t %f",
1015  iDet,
1016  parameters[iDet*fgNParCh+0],parameters[iDet*fgNParCh+1],
1017  parameters[iDet*fgNParCh+3],parameters[iDet*fgNParCh+2]
1018  ) );
1019  }
1020 
1021 }
1022 
1023 //_____________________________________________________
1026 
1027 //_____________________________________________________
1028 Double_t AliMUONAlignment::GetParError(Int_t iPar) const
1029 { return fMillepede->GetParError(iPar); }
1030 
1031 //______________________________________________________________________
1033  const AliMUONGeometryTransformer * transformer,
1034  const double *misAlignments, Bool_t )
1035 {
1036 
1039 
1040  // Takes the internal geometry module transformers, copies them
1041  // and gets the Detection Elements from them.
1042  // Takes misalignment parameters and applies these
1043  // to the local transform of the Detection Element
1044  // Obtains the global transform by multiplying the module transformer
1045  // transformation with the local transformation
1046  // Applies the global transform to a new detection element
1047  // Adds the new detection element to a new module transformer
1048  // Adds the new module transformer to a new geometry transformer
1049  // Returns the new geometry transformer
1050 
1051  Double_t lModuleMisAlignment[fgNParCh] = {0};
1052  Double_t lDetElemMisAlignment[fgNParCh] = {0};
1053  const TClonesArray* oldMisAlignArray( transformer->GetMisAlignmentData() );
1054 
1055  AliMUONGeometryTransformer *newGeometryTransformer = new AliMUONGeometryTransformer();
1056  for( Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); ++iMt )
1057  {
1058 
1059  // module transformers
1060  const AliMUONGeometryModuleTransformer *kModuleTransformer = transformer->GetModuleTransformer(iMt, kTRUE);
1061 
1062  AliMUONGeometryModuleTransformer *newModuleTransformer = new AliMUONGeometryModuleTransformer(iMt);
1063  newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1064 
1065  // get transformation
1066  TGeoHMatrix deltaModuleTransform( DeltaTransform( lModuleMisAlignment ) );
1067 
1068  // update module
1069  TGeoHMatrix moduleTransform( *kModuleTransformer->GetTransformation() );
1070  TGeoHMatrix newModuleTransform( AliMUONGeometryBuilder::Multiply( deltaModuleTransform, moduleTransform ) );
1071  newModuleTransformer->SetTransformation(newModuleTransform);
1072 
1073  // Get matching old alignment and update current matrix accordingly
1074  if( oldMisAlignArray )
1075  {
1076 
1077  const AliAlignObjMatrix* oldAlignObj(0);
1078  const Int_t moduleId( kModuleTransformer->GetModuleId() );
1079  const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, moduleId );
1080  for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
1081  {
1082 
1083  const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
1084  if( localAlignObj && localAlignObj->GetVolUID() == volId )
1085  {
1086  oldAlignObj = localAlignObj;
1087  break;
1088  }
1089 
1090  }
1091 
1092  // multiply
1093  if( oldAlignObj )
1094  {
1095 
1096  TGeoHMatrix oldMatrix;
1097  oldAlignObj->GetMatrix( oldMatrix );
1098  deltaModuleTransform.Multiply( &oldMatrix );
1099 
1100  }
1101 
1102  }
1103 
1104  // Create module mis alignment matrix
1105  newGeometryTransformer->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
1106 
1107  AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
1108 
1109  TIter next(detElements->CreateIterator());
1110  AliMUONGeometryDetElement* detElement;
1111  Int_t iDe(-1);
1112  while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
1113  {
1114  ++iDe;
1115  // make a new detection element
1116  AliMUONGeometryDetElement *newDetElement = new AliMUONGeometryDetElement(detElement->GetId(), detElement->GetVolumePath());
1117  TString lDetElemName(detElement->GetDEName());
1118  lDetElemName.ReplaceAll("DE","");
1119 
1120  // store detector element id and number
1121  const Int_t iDetElemId = lDetElemName.Atoi();
1122  if( DetElemIsValid( iDetElemId ) )
1123  {
1124 
1125  const Int_t iDetElemNumber( GetDetElemNumber( iDetElemId ) );
1126 
1127  for( int i=0; i<fgNParCh; ++i )
1128  {
1129  lDetElemMisAlignment[i] = 0.0;
1130  if( iMt<fgNTrkMod ) { lDetElemMisAlignment[i] = misAlignments[iDetElemNumber*fgNParCh+i]; }
1131  }
1132 
1133  // get transformation
1134  TGeoHMatrix deltaGlobalTransform( DeltaTransform( lDetElemMisAlignment ) );
1135 
1136  // update module
1137  TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
1138  TGeoHMatrix newGlobalTransform( AliMUONGeometryBuilder::Multiply( deltaGlobalTransform, globalTransform ) );
1139  newDetElement->SetGlobalTransformation( newGlobalTransform );
1140  newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
1141 
1142  // Get matching old alignment and update current matrix accordingly
1143  if( oldMisAlignArray )
1144  {
1145 
1146  const AliAlignObjMatrix* oldAlignObj(0);
1147  const int detElemId( detElement->GetId() );
1148  const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
1149  for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
1150  {
1151 
1152  const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
1153  if( localAlignObj && localAlignObj->GetVolUID() == volId )
1154  {
1155  oldAlignObj = localAlignObj;
1156  break;
1157  }
1158 
1159  }
1160 
1161  // multiply
1162  if( oldAlignObj )
1163  {
1164 
1165  TGeoHMatrix oldMatrix;
1166  oldAlignObj->GetMatrix( oldMatrix );
1167  deltaGlobalTransform.Multiply( &oldMatrix );
1168 
1169  }
1170 
1171  }
1172 
1173  // Create misalignment matrix
1174  newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
1175 
1176  } else {
1177 
1178  // "invalid" detector elements come from MTR and are left unchanged
1179  AliInfo( Form( "Keeping detElement %i unchanged", iDetElemId ) );
1180 
1181  // update module
1182  TGeoHMatrix globalTransform( *detElement->GetGlobalTransformation() );
1183  newDetElement->SetGlobalTransformation( globalTransform );
1184  newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(), newDetElement);
1185 
1186  // Get matching old alignment and update current matrix accordingly
1187  if( oldMisAlignArray )
1188  {
1189 
1190  const AliAlignObjMatrix* oldAlignObj(0);
1191  const int detElemId( detElement->GetId() );
1192  const Int_t volId = AliGeomManager::LayerToVolUID(AliGeomManager::kMUON, detElemId );
1193  for( Int_t pos = 0; pos < oldMisAlignArray->GetEntriesFast(); ++pos )
1194  {
1195 
1196  const AliAlignObjMatrix* localAlignObj( dynamic_cast<const AliAlignObjMatrix*>(oldMisAlignArray->At( pos ) ) );
1197  if( localAlignObj && localAlignObj->GetVolUID() == volId )
1198  {
1199  oldAlignObj = localAlignObj;
1200  break;
1201  }
1202 
1203  }
1204 
1205  // multiply
1206  if( oldAlignObj )
1207  {
1208 
1209  TGeoHMatrix oldMatrix;
1210  oldAlignObj->GetMatrix( oldMatrix );
1211  newGeometryTransformer->AddMisAlignDetElement(detElement->GetId(), oldMatrix );
1212 
1213  }
1214 
1215  }
1216 
1217  }
1218 
1219  }
1220 
1221  newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1222 
1223  }
1224 
1225  return newGeometryTransformer;
1226 
1227 }
1228 
1229 //______________________________________________________________________
1230 void AliMUONAlignment::SetAlignmentResolution( const TClonesArray* misAlignArray, Int_t rChId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY )
1231 {
1232 
1235  TMatrixDSym mChCorrMatrix(6);
1236  mChCorrMatrix[0][0]=chResX*chResX;
1237  mChCorrMatrix[1][1]=chResY*chResY;
1238 
1239  TMatrixDSym mDECorrMatrix(6);
1240  mDECorrMatrix[0][0]=deResX*deResX;
1241  mDECorrMatrix[1][1]=deResY*deResY;
1242 
1243  AliAlignObjMatrix *alignMat = 0x0;
1244 
1245  for( Int_t chId = 0; chId <= 9; ++chId )
1246  {
1247 
1248  // skip chamber if selection is valid, and does not match
1249  if( rChId > 0 && chId+1 != rChId ) continue;
1250 
1251  TString chName1;
1252  TString chName2;
1253  if (chId<4)
1254  {
1255 
1256  chName1 = Form("GM%d",chId);
1257  chName2 = Form("GM%d",chId);
1258 
1259  } else {
1260 
1261  chName1 = Form("GM%d",4+(chId-4)*2);
1262  chName2 = Form("GM%d",4+(chId-4)*2+1);
1263 
1264  }
1265 
1266  for( int i=0; i<misAlignArray->GetEntries(); ++i )
1267  {
1268 
1269  alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
1270  TString volName(alignMat->GetSymName());
1271  if((volName.Contains(chName1)&&
1272  ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
1273  (volName.Length()==volName.Index(chName1)+chName1.Length())))||
1274  (volName.Contains(chName2)&&
1275  ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
1276  (volName.Length()==volName.Index(chName2)+chName2.Length()))))
1277  {
1278 
1279  volName.Remove(0,volName.Last('/')+1);
1280  if (volName.Contains("GM")) alignMat->SetCorrMatrix(mChCorrMatrix);
1281  else if (volName.Contains("DE")) alignMat->SetCorrMatrix(mDECorrMatrix);
1282 
1283  }
1284 
1285  }
1286 
1287  }
1288 
1289 }
1290 
1291 //_____________________________________________________
1293 {
1294 
1295  // initialize matrices
1296  TMatrixD AtGASum(4,4);
1297  AtGASum.Zero();
1298 
1299  TMatrixD AtGMSum(4,1);
1300  AtGMSum.Zero();
1301 
1302  // loop over clusters
1303  TObjArray* trackParamAtCluster( track->GetTrackParamAtCluster() );
1304  for( Int_t iTrackParam = 0; iTrackParam < trackParamAtCluster->GetEntries(); ++iTrackParam )
1305  {
1306 
1307  // get track parameters
1308  AliMUONTrackParam* trackParam( (AliMUONTrackParam *) trackParamAtCluster->At(iTrackParam) );
1309  if( !trackParam ) continue;
1310 
1311  // get cluster
1312  AliMUONVCluster* cluster = trackParam->GetClusterPtr();
1313  if( !cluster ) continue;
1314 
1315  // projection matrix
1316  TMatrixD A(2, 4);
1317  A.Zero();
1318  A(0,0) = 1; A(0, 2) = (cluster->GetZ() - z0);
1319  A(1,1) = 1; A(1, 3) = (cluster->GetZ() - z0);
1320 
1321  TMatrixD At( TMatrixD::kTransposed, A );
1322 
1323  // gain matrix
1324  TMatrixD G(2,2);
1325  G.Zero();
1326  G(0,0) = 1.0/Square(cluster->GetErrX());
1327  G(1,1) = 1.0/Square(cluster->GetErrY());
1328 
1329  const TMatrixD AtG( At, TMatrixD::kMult, G );
1330  const TMatrixD AtGA( AtG, TMatrixD::kMult, A );
1331  AtGASum += AtGA;
1332 
1333  // measurement
1334  TMatrixD M(2,1);
1335  M(0,0) = cluster->GetX();
1336  M(1,0) = cluster->GetY();
1337  const TMatrixD AtGM( AtG, TMatrixD::kMult, M );
1338  AtGMSum += AtGM;
1339 
1340  }
1341 
1342  // perform inversion
1343  TMatrixD AtGASumInv( TMatrixD::kInverted, AtGASum );
1344  TMatrixD X( AtGASumInv, TMatrixD::kMult, AtGMSum );
1345 
1346 // // TODO: compare with initial track parameters
1347 // AliInfo( Form( "x: %.3f vs %.3f", fTrackPos0[0], X(0,0) ) );
1348 // AliInfo( Form( "y: %.3f vs %.3f", fTrackPos0[1], X(1,0) ) );
1349 // AliInfo( Form( "dxdz: %.6g vs %.6g", fTrackSlope0[0], X(2,0) ) );
1350 // AliInfo( Form( "dydz: %.6g vs %.6g\n", fTrackSlope0[1], X(3,0) ) );
1351 
1352  // fill output parameters
1353  LocalTrackParam out;
1354  out.fTrackX = X(0,0);
1355  out.fTrackY = X(1,0);
1356  out.fTrackZ = z0;
1357  out.fTrackSlopeX = X(2,0);
1358  out.fTrackSlopeY = X(3,0);
1359 
1360  return out;
1361 
1362 }
1363 
1364 //_____________________________________________________
1366 {
1367 
1369  // get detector element number from Alice ID
1370  const Int_t detElemId = cluster->GetDetElemId();
1371  fDetElemNumber = GetDetElemNumber( detElemId );
1372 
1373  // get detector element
1374  const AliMUONGeometryDetElement* detElement = fTransform->GetDetElement( detElemId );
1375 
1376  /*
1377  get the global transformation matrix and store its inverse, in order to manually perform
1378  the global to Local transformations needed to calculate the derivatives
1379  */
1380  fGeoCombiTransInverse = detElement->GetGlobalTransformation()->Inverse();
1381 
1382 }
1383 
1384 //______________________________________________________________________
1386 {
1387 
1389  fClustPos[0] = cluster->GetX();
1390  fClustPos[1] = cluster->GetY();
1391  fClustPos[2] = cluster->GetZ();
1392 
1393 }
1394 
1395 //______________________________________________________________________
1397 {
1398 
1400  fTrackPos[0] = trackParam->GetNonBendingCoor();
1401  fTrackPos[1] = trackParam->GetBendingCoor();
1402  fTrackPos[2] = trackParam->GetZ();
1403  fTrackSlope[0] = trackParam->GetNonBendingSlope();
1404  fTrackSlope[1] = trackParam->GetBendingSlope();
1405 
1406 }
1407 
1408 //______________________________________________________________________
1410 {
1412 
1413  // 'inverse' (GlobalToLocal) rotation matrix
1414  const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
1415 
1416  // local derivatives
1417  SetLocalDerivative( 0, r[0] );
1418  SetLocalDerivative( 1, r[0]*(fTrackPos[2] - fTrackPos0[2]) );
1419  SetLocalDerivative( 2, r[1] );
1420  SetLocalDerivative( 3, r[1]*(fTrackPos[2] - fTrackPos0[2]) );
1421 
1422  // global derivatives
1423  /*
1424  alignment parameters are
1425  0: delta_x
1426  1: delta_y
1427  2: delta_phiz
1428  3: delta_z
1429  */
1430 
1433 
1434  if( fBFieldOn )
1435  {
1436 
1437  // use local position for derivatives vs 'delta_phi_z'
1438  SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*fTrackPos[0] + r[0]*fTrackPos[1] );
1439 
1440  // use local slopes for derivatives vs 'delta_z'
1442 
1443  } else {
1444 
1445  // local copy of extrapolated track positions
1446  const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1447  const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
1448 
1449  // use properly extrapolated position for derivatives vs 'delta_phi_z'
1450  SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[1]*trackPosX + r[0]*trackPosY );
1451 
1452  // use slopes at origin for derivatives vs 'delta_z'
1454 
1455  }
1456 
1457  // store local equation
1459 
1460 }
1461 
1462 //______________________________________________________________________
1464 {
1466 
1467  // 'inverse' (GlobalToLocal) rotation matrix
1468  const Double_t* r( fGeoCombiTransInverse.GetRotationMatrix() );
1469 
1470  // store local derivatives
1471  SetLocalDerivative( 0, r[3] );
1472  SetLocalDerivative( 1, r[3]*(fTrackPos[2] - fTrackPos0[2] ) );
1473  SetLocalDerivative( 2, r[4] );
1474  SetLocalDerivative( 3, r[4]*(fTrackPos[2] - fTrackPos0[2] ) );
1475 
1476  // set global derivatives
1479 
1480  if( fBFieldOn )
1481  {
1482 
1483  // use local position for derivatives vs 'delta_phi'
1485 
1486  // use local slopes for derivatives vs 'delta_z'
1488 
1489  } else {
1490 
1491  // local copy of extrapolated track positions
1492  const Double_t trackPosX = fTrackPos0[0]+fTrackSlope0[0]*( fTrackPos[2]-fTrackPos0[2] );
1493  const Double_t trackPosY = fTrackPos0[1]+fTrackSlope0[1]*( fTrackPos[2]-fTrackPos0[2] );
1494 
1495  // use properly extrapolated position for derivatives vs 'delta_phi'
1496  SetGlobalDerivative( fDetElemNumber*fgNParCh + 2, -r[4]*trackPosX + r[3]*trackPosY );
1497 
1498  // use slopes at origin for derivatives vs 'delta_z'
1500 
1501  }
1502 
1503  // store local equation
1505 
1506 }
1507 
1508 //_________________________________________________________________________
1509 TGeoCombiTrans AliMUONAlignment::DeltaTransform( const double *lMisAlignment) const
1510 {
1512 
1513  // translation
1514  const TGeoTranslation deltaTrans( lMisAlignment[0], lMisAlignment[1], lMisAlignment[3]);
1515 
1516  // rotation
1517  TGeoRotation deltaRot;
1518  deltaRot.RotateZ(lMisAlignment[2]*180./TMath::Pi());
1519 
1520  // combined rotation and translation.
1521  return TGeoCombiTrans(deltaTrans,deltaRot);
1522 
1523 }
1524 
1525 //______________________________________________________________________
1526 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value)
1527 {
1529  if( !fInitialized )
1530  { AliFatal( "Millepede is not initialized" ); }
1531 
1532  fMillepede->SetGlobalConstraint(par, value);
1533 }
1534 
1535 //______________________________________________________________________
1536 Bool_t AliMUONAlignment::DetElemIsValid( Int_t iDetElemId ) const
1537 {
1539  const Int_t iCh = iDetElemId/100;
1540  const Int_t iDet = iDetElemId%100;
1541  return ( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] );
1542 }
1543 
1544 //______________________________________________________________________
1545 Int_t AliMUONAlignment::GetDetElemNumber( Int_t iDetElemId ) const
1546 {
1548  // get chamber and element number in chamber
1549  const Int_t iCh = iDetElemId/100;
1550  const Int_t iDet = iDetElemId%100;
1551 
1552  // make sure detector index is valid
1553  if( !( iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh-1] ) )
1554  { AliFatal( Form( "Invalid detector element id: %i", iDetElemId ) ); }
1555 
1556  // add number of detectors up to this chamber
1557  return iDet + fgSNDetElemCh[iCh-1];
1558 
1559 }
1560 
1561 //______________________________________________________________________
1562 Int_t AliMUONAlignment::GetChamberId( Int_t iDetElemNumber ) const
1563 {
1565  Int_t iCh( 0 );
1566  for( iCh=0; iCh<fgNCh; iCh++ )
1567  { if( iDetElemNumber < fgSNDetElemCh[iCh] ) break; }
1568 
1569  return iCh;
1570 }
1571 
1572 //______________________________________________________________________
1573 TString AliMUONAlignment::GetParameterMaskString( UInt_t mask ) const
1574 {
1575  TString out;
1576  if( mask & ParX ) out += "X";
1577  if( mask & ParY ) out += "Y";
1578  if( mask & ParZ ) out += "Z";
1579  if( mask & ParTZ ) out += "T";
1580  return out;
1581 }
1582 
1583 //______________________________________________________________________
1584 TString AliMUONAlignment::GetSidesMaskString( UInt_t mask ) const
1585 {
1586  TString out;
1587  if( mask & SideTop ) out += "T";
1588  if( mask & SideLeft ) out += "L";
1589  if( mask & SideBottom ) out += "B";
1590  if( mask & SideRight ) out += "R";
1591  return out;
1592 }
void GroupChamber(Int_t iCh, UInt_t parameterMask=ParAll)
Class for alignment of muon spectrometer.
void InitGlobalParameters(Double_t *par)
initialize global parameters to a give set of values
Array(void)
contructor
self initialized array, used for adding constraints
Array & operator=(const Array &)
Not implemented.
const AliMUONGeometryDetElement * GetDetElement(Int_t detElemId, Bool_t warn=true) const
void SetCorrMatrix(Double_t *cov)
Number tracking chambers.
Int_t GetNofModuleTransformers() const
Return the number of contained module transformers.
Geometry transformer for a detector module.
virtual Double_t GetZ() const =0
Return coordinate Z (cm)
#define TObjArray
Int_t fNStdDev
Number of standard deviations for chi2 cut.
void SetChamberNonLinear(Int_t iCh, UInt_t parameterMask)
Double_t GetBendingCoor() const
return bending coordinate (cm)
void FixAll(UInt_t parameterMask=ParAll)
Top container class for geometry transformations.
Int_t fGlobalParameterStatus[fNGlobal]
void SetGlobalDerivative(Int_t index, Double_t value)
Set array of global derivatives.
void SetRecordRun(Int_t run)
void FixParameter(Int_t iPar)
Double_t fLocalDerivatives[fNLocal]
Array of local derivatives.
LocalTrackParam * fTrackParamNew
void SetGlobalParameters(Double_t *par)
virtual Double_t GetErrX() const =0
Return resolution (cm) on coordinate X.
Double_t GetZ() const
return Z coordinate (cm)
Double_t fTrackSlope0[2]
Track slope at reference point.
Double_t fTrackSlope[2]
Track slope at current point.
Track parameters in ALICE dimuon spectrometer.
void SetGlobalConstraint(const double *dergb, double val, double sigma=0)
static int fgNDetElemCh[10]
Int_t GetId() const
Return detection element ID.
AliMUONGeometryTransformer * fTransform
Geometry transformation.
Bool_t InitDataRecStorage(Bool_t read=kFALSE)
TFile * fTFile
output TFile
AliTPCfastTrack * track
LocalTrackParam RefitStraightTrack(AliMUONTrack *, Double_t) const
refit track using straight track model
void GroupHalfChamber(Int_t iCh, Int_t iHalf, UInt_t parameterMask=ParAll)
static const Int_t fgNDetElemCh[fgNCh]
Number of detection elements per chamber.
Int_t fDetElemNumber
current detection element number
TTree * fTTree
output TTree
void AddConstraints(const Bool_t *bChOnOff, UInt_t parameterMask)
virtual void GetMatrix(TGeoHMatrix &m) const
void SetNonLinear(int index, Bool_t v=kTRUE)
void SetRecordWeight(double wgh)
void ReleaseChamber(Int_t iCh, UInt_t parameterMask=ParAll)
void FixHalfSpectrometer(const Bool_t *bChOnOff, UInt_t sidesMask=AllSides, UInt_t parameterMask=ParAll)
#define M(PID)
Definition: AliPID.cxx:49
void ReleaseDetElem(Int_t iDetElemId, UInt_t parameterMask=ParAll)
TString GetSidesMaskString(UInt_t sidesMask) const
void SaveRecordData()
const TClonesArray * GetMisAlignmentData() const
Return the array of misalignment data.
void ReleaseParameter(Int_t iPar)
void SetGlobalTransformation(const TGeoHMatrix &transform, Bool_t warn=true)
Double_t GetBendingSlope() const
return bending slope (cm ** -1)
AliMillePedeRecord * GetRecord() const
void Local2Global(Int_t detElemId, Float_t xl, Float_t yl, Float_t zl, Float_t &xg, Float_t &yg, Float_t &zg) const
void SetParameterNonLinear(Int_t iPar)
AliMillePede2 * fMillepede
Detector independent alignment class.
void PrintGlobalParameters(void) const
print global parameters
abstract base class for clusters
void AddConstraint(Double_t *parameters, Double_t value)
void GlobalFit(Double_t *parameters, Double_t *errors, Double_t *pulls)
perform global fit
Double_t GetParError(int iPar) const
local track parameters, for refit
max number of detector elements per half chamber
TGeoCombiTrans DeltaTransform(const double *detElemMisAlignment) const
Number of degrees of freedom per chamber.
Number of local parameters.
Number of half chambers.
#define AliInfo(message)
Definition: AliLog.h:484
const TGeoHMatrix * GetTransformation() const
Return the module transformation wrt to the top volume (world)
Double_t fStartFac
Initial value for chi2 cut.
AliMillePedeRecord * ProcessTrack(AliMUONTrack *track, Bool_t doAlignment, Double_t weight=1)
Bool_t fBFieldOn
Flag for Magnetic filed On/Off.
virtual Double_t GetErrY() const =0
Return resolution (cm) on coordinate Y.
Int_t InitMille(int nGlo, int nLoc, Int_t lNStdDev=-1, double lResCut=-1., double lResCutInit=-1., const Int_t *regroup=0)
AliMillePedeRecord fTrackRecord
running Track record
static TGeoHMatrix Multiply(const TGeoMatrix &m1, const TGeoMatrix &m2)
Int_t GetChamberId(Int_t iDetElemNumber) const
const TGeoHMatrix * GetGlobalTransformation() const
Return the detection element transformation wrt world.
Bool_t fDoEvaluation
preform evaluation
Number of global parameters.
Double_t fResCut
Cut on residual for other iterations.
Double_t GetParError(Int_t iPar) const
get error on a given parameter
AliMUONGeometryTransformer * ReAlign(const AliMUONGeometryTransformer *transformer, const double *misAlignments, Bool_t verbose)
void GroupDetElems(Int_t detElemMin, Int_t detElemMax, UInt_t parameterMask=ParAll)
Int_t fRunNumber
current run id
void Add(Int_t keyFirst, Int_t keySecond, TObject *object)
Definition: AliMpExMap.cxx:292
void SetParSigma(Int_t i, Double_t par)
static const Int_t fgNDetElemHalfCh[fgNHalfCh]
Number of detection element per tracking module.
Number of tracking modules.
#define AliFatal(message)
Definition: AliLog.h:640
Double_t fTrackPos[3]
Track intersection at current point.
void FillRecPointData(AliMUONVCluster *)
TString GetParameterMaskString(UInt_t parameterMask) const
Int_t GetModuleId() const
Return module ID.
TVectorD errors
Definition: driftITSTPC.C:97
void SetDetElemNonLinear(Int_t iSt, UInt_t parameterMask)
virtual Double_t GetY() const =0
Return coordinate Y (cm)
static const Int_t fgDetElemHalfCh[fgNHalfCh][fgNDetHalfChMax]
list of detection elements per tracking module
void SetSigmaXY(Double_t sigmaX, Double_t sigmaY)
AliMpExMap * GetDetElementStore() const
Return detection elements associated with this module.
void SetAlignmentResolution(const TClonesArray *misAlignArray, Int_t chId, Double_t chResX, Double_t chResY, Double_t deResX, Double_t deResY)
Double_t fGlobalDerivatives[fNGlobal]
Array of global derivatives.
void FixDetElem(Int_t iDetElemId, UInt_t parameterMask=ParAll)
void FixChamber(Int_t iCh, UInt_t parameterMask=ParAll)
Double_t fMeas[2]
Current measurement (depend on B field On/Off)
Double_t GetNonBendingCoor() const
return non bending coordinate (cm)
Bool_t fInitialized
true when initialized
Double_t Square(Double_t x)
static UShort_t LayerToVolUID(ELayerID layerId, Int_t modId)
Double_t values[AliMUONAlignment::fNGlobal]
array
static Int_t GetDetElemId(UInt_t uniqueID)
Return detection element id, part of the uniqueID.
const char * GetSymName() const
Definition: AliAlignObj.h:63
void FillDetElemData(AliMUONVCluster *)
UShort_t GetVolUID() const
Definition: AliAlignObj.h:64
virtual Double_t GetX() const =0
Return coordinate X (cm)
Double_t fAllowVar[fgNParCh]
"Encouraged" variation for degrees of freedom
void SetLocalEquation(double *dergb, double *derlc, double lMeas, double lSigma)
void FillTrackParamData(AliMUONTrackParam *)
Int_t SetIterations(double lChi2CutFac)
const AliMUONGeometryModuleTransformer * GetModuleTransformer(Int_t index, Bool_t warn=true) const
Reconstructed track in ALICE dimuon spectrometer.
Definition: AliMUONTrack.h:24
void SetAllowedVariation(Int_t iPar, Double_t value)
Int_t GetDetElemNumber(Int_t iDetElemId) const
Double_t fTrackPos0[3]
Track intersection at reference point.
Bool_t fRefitStraightTracks
true if straight track refit is to be performed
void SetTransformation(const TGeoHMatrix &transform)
Bool_t DetElemIsValid(Int_t iDetElemId) const
Class for storing detection element transformations.
Double_t fSigma[2]
Estimated resolution on measurement.
LocalTrackParam * fTrackParamOrig
original local track params
TGeoCombiTrans fGeoCombiTransInverse
Int_t GlobalFit(Double_t *par=0, Double_t *error=0, Double_t *pull=0)
Double_t fClustPos[3]
Cluster (global) position.
class TMatrixT< Double_t > TMatrixD
Helper class making Root persistent TExMap.
Definition: AliMpExMap.h:28
AliMpExMapIterator * CreateIterator() const
Definition: AliMpExMap.cxx:357
void SetLocalDerivative(Int_t index, Double_t value)
Set array of local derivatives.
Double_t fResCutInitial
Cut on residual for first iteration.
static const Int_t fgSNDetElemCh[fgNCh+1]
Sum of detection elements up to this chamber.
TObjArray * GetTrackParamAtCluster() const
Int_t PrintGlobalParameters() const
Double_t GetNonBendingSlope() const
return non bending slope (cm ** -1)