1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONAlignmentV5
20 /// Alignment class fro the ALICE DiMuon spectrometer
22 /// MUON specific alignment class which interface to AliMillepede.
23 /// For each track ProcessTrack calculates the local and global derivatives
24 /// at each hit and fill the corresponding local equations. Provide methods for
25 /// fixing or constraining detection elements for best results.
27 /// \author Bruce Becker, Javier Castillo
28 //-----------------------------------------------------------------------------
30 #include "AliMUONAlignment.h"
31 #include "AliMUONTrack.h"
32 #include "AliMUONRawCluster.h"
33 #include "AliMUONTrackParam.h"
34 #include "AliMUONHitForRec.h"
35 #include "AliMUONGeometryTransformer.h"
36 #include "AliMUONGeometryModuleTransformer.h"
37 #include "AliMUONGeometryDetElement.h"
38 #include "AliMUONGeometryStore.h"
39 #include "AliMUONGeometryBuilder.h"
42 #include "AliMUONConstants.h"
44 #include "AliMillepede.h"
46 ClassImp(AliMUONAlignment)
47 Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
48 Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
49 Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
50 Int_t AliMUONAlignment::fgNParCh = 3;
51 Int_t AliMUONAlignment::fgNCh = 10;
52 Int_t AliMUONAlignment::fgNSt = 5;
54 AliMUONAlignment::AliMUONAlignment()
66 fNGlobal(fgNDetElem*fgNParCh),
76 /// Default constructor, setting define alignment parameters
80 fDoF[0] = kTRUE; fDoF[1] = kTRUE; fDoF[2] = kTRUE;
81 fAllowVar[0] = 0.05; fAllowVar[1] = 0.05; fAllowVar[2] = 0.001;
83 AliInfo(Form("fAllowVar[0]: %f\t fAllowVar[1]: %f\t fPhi: %f\t fgNDetElem: %i\t fNGlobal: %i\t fNLocal: %i",fAllowVar[0],fAllowVar[1],fPhi,fgNDetElem,fNGlobal,fNLocal));
85 fMillepede = new AliMillepede();
87 Init(fNGlobal, fNLocal, fNStdDev);
90 AliInfo("Parameters initialized to zero");
94 AliMUONAlignment::~AliMUONAlignment() {
98 void AliMUONAlignment::Init(Int_t nGlobal, /* number of global paramers */
99 Int_t nLocal, /* number of local parameters */
100 Int_t nStdDev /* std dev cut */ )
102 /// Initialization of AliMillepede. Fix parameters, define constraints ...
103 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
105 Bool_t bStOnOff[5] = {kFALSE,kFALSE,kTRUE,kTRUE,kTRUE};
107 AllowVariations(bStOnOff);
109 // Fix parameters or add constraints here
110 for (Int_t iSt=0; iSt<5; iSt++)
111 if (!bStOnOff[iSt]) FixStation(iSt+1);
113 Bool_t bVarXYT[3] = {kFALSE,kTRUE,kFALSE};
114 Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
116 AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
117 bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
118 bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
119 AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
120 bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
121 AddConstraints(bStOnOff,bVarXYT);
124 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
127 void AliMUONAlignment::FixStation(Int_t iSt){
128 /// Fix all detection elements of station iSt
129 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
130 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
131 for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){
132 FixParameter(i*fgNParCh+0, 0.0);
133 FixParameter(i*fgNParCh+1, 0.0);
134 FixParameter(i*fgNParCh+2, 0.0);
138 void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){
139 /// Set non linear parameter flag selected stations and degrees of freedom
140 for (Int_t i = 0; i < fgNDetElem; i++){
142 for (iCh=1; iCh<=fgNCh; iCh++){
143 if (i<fgSNDetElemCh[iCh-1]) break;
145 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
147 if (lVarXYT[0]) { // X constraint
148 SetNonLinear(i*fgNParCh+0);
150 if (lVarXYT[1]) { // Y constraint
151 SetNonLinear(i*fgNParCh+1);
153 if (lVarXYT[2]) { // T constraint
154 SetNonLinear(i*fgNParCh+2);
160 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT){
161 /// Add constraint equations for selected stations and degrees of freedom
162 for (Int_t i = 0; i < fgNDetElem; i++){
164 for (iCh=1; iCh<=fgNCh; iCh++){
165 if (i<fgSNDetElemCh[iCh-1]) break;
167 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
169 if (lVarXYT[0]) { // X constraint
170 fConstraintX[i*fgNParCh+0]=1.0;
172 if (lVarXYT[1]) { // Y constraint
173 fConstraintY[i*fgNParCh+1]=1.0;
175 if (lVarXYT[2]) { // T constraint
176 fConstraintP[i*fgNParCh+2]=1.0;
180 if (lVarXYT[0]) { // X constraint
181 AddConstraint(fConstraintX,0.0);
183 if (lVarXYT[1]) { // Y constraint
184 AddConstraint(fConstraintY,0.0);
186 if (lVarXYT[2]) { // T constraint
187 AddConstraint(fConstraintP,0.0);
191 void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t *lDetTLBR){
192 /// Add constraint equations for selected stations, degrees of freedom detector half
193 for (Int_t i = 0; i < fgNDetElem; i++){
195 for (iCh=1; iCh<=fgNCh; iCh++){
196 if (i<fgSNDetElemCh[iCh-1]) break;
198 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
200 if (lVarXYT[0]) { // X constraint
201 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
202 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
203 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
204 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
206 if (lVarXYT[1]) { // X constraint
207 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
208 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
209 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
210 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
212 if (lVarXYT[2]) { // X constraint
213 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
214 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
215 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
216 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
220 if (lVarXYT[0]) { // X constraint
221 if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
222 if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
223 if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
224 if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
226 if (lVarXYT[1]) { // Y constraint
227 if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
228 if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
229 if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
230 if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
232 if (lVarXYT[2]) { // T constraint
233 if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
234 if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
235 if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
236 if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
240 void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar){
241 /// Set constrain equation for top half of spectrometer
242 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
243 if (lCh>=1 && lCh<=4){
244 if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
245 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
248 if (lCh>=5 && lCh<=6){
249 if (lDetElemNumber>=0&&lDetElemNumber<=9){
250 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
253 if (lCh>=7 && lCh<=10){
254 if (lDetElemNumber>=0&&lDetElemNumber<=13){
255 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
260 void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar){
261 /// Set constrain equation for left half of spectrometer
262 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
263 if (lCh>=1 && lCh<=4){
264 if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
265 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
268 if (lCh>=5 && lCh<=6){
269 if (lDetElemNumber>=5&&lDetElemNumber<=13){
270 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
273 if (lCh>=7 && lCh<=10){
274 if (lDetElemNumber>=7&&lDetElemNumber<=19){
275 lConstraintL[lDetElem*fgNParCh+iVar]=1.0;
280 void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar){
281 /// Set constrain equation for bottom half of spectrometer
282 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
283 if (lCh>=1 && lCh<=4){
284 if (lDetElemNumber==2 && lDetElemNumber==3){ // From track crossings
285 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
288 if (lCh>=5 && lCh<=6){
289 if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
290 (lDetElemNumber==0)){
291 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
294 if (lCh>=7 && lCh<=10){
295 if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
296 (lDetElemNumber==0)){
297 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
302 void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar){
303 /// Set constrain equation for right half of spectrometer
304 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
305 if (lCh>=1 && lCh<=4){
306 if (lDetElemNumber==0 && lDetElemNumber==3){ // From track crossings
307 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
310 if (lCh>=5 && lCh<=6){
311 if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
312 (lDetElemNumber>=14&&lDetElemNumber<=17)){
313 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
316 if (lCh>=7 && lCh<=10){
317 if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
318 (lDetElemNumber>=20&&lDetElemNumber<=25)){
319 lConstraintR[lDetElem*fgNParCh+iVar]=1.0;
324 void AliMUONAlignment::ResetConstraints(){
325 /// Reset all constraint equations
326 for (Int_t i = 0; i < fgNDetElem; i++){
327 fConstraintX[i*fgNParCh+0]=0.0;
328 fConstraintX[i*fgNParCh+1]=0.0;
329 fConstraintX[i*fgNParCh+2]=0.0;
330 fConstraintY[i*fgNParCh+0]=0.0;
331 fConstraintY[i*fgNParCh+1]=0.0;
332 fConstraintY[i*fgNParCh+2]=0.0;
333 fConstraintP[i*fgNParCh+0]=0.0;
334 fConstraintP[i*fgNParCh+1]=0.0;
335 fConstraintP[i*fgNParCh+2]=0.0;
336 fConstraintXT[i*fgNParCh+0]=0.0;
337 fConstraintXT[i*fgNParCh+1]=0.0;
338 fConstraintXT[i*fgNParCh+2]=0.0;
339 fConstraintYT[i*fgNParCh+0]=0.0;
340 fConstraintYT[i*fgNParCh+1]=0.0;
341 fConstraintYT[i*fgNParCh+2]=0.0;
342 fConstraintPT[i*fgNParCh+0]=0.0;
343 fConstraintPT[i*fgNParCh+1]=0.0;
344 fConstraintPT[i*fgNParCh+2]=0.0;
345 fConstraintXL[i*fgNParCh+0]=0.0;
346 fConstraintXL[i*fgNParCh+1]=0.0;
347 fConstraintXL[i*fgNParCh+2]=0.0;
348 fConstraintYL[i*fgNParCh+0]=0.0;
349 fConstraintYL[i*fgNParCh+1]=0.0;
350 fConstraintYL[i*fgNParCh+2]=0.0;
351 fConstraintPL[i*fgNParCh+0]=0.0;
352 fConstraintPL[i*fgNParCh+1]=0.0;
353 fConstraintPL[i*fgNParCh+2]=0.0;
354 fConstraintXB[i*fgNParCh+0]=0.0;
355 fConstraintXB[i*fgNParCh+1]=0.0;
356 fConstraintXB[i*fgNParCh+2]=0.0;
357 fConstraintYB[i*fgNParCh+0]=0.0;
358 fConstraintYB[i*fgNParCh+1]=0.0;
359 fConstraintYB[i*fgNParCh+2]=0.0;
360 fConstraintPB[i*fgNParCh+0]=0.0;
361 fConstraintPB[i*fgNParCh+1]=0.0;
362 fConstraintPB[i*fgNParCh+2]=0.0;
363 fConstraintXR[i*fgNParCh+0]=0.0;
364 fConstraintXR[i*fgNParCh+1]=0.0;
365 fConstraintXR[i*fgNParCh+2]=0.0;
366 fConstraintYR[i*fgNParCh+0]=0.0;
367 fConstraintYR[i*fgNParCh+1]=0.0;
368 fConstraintYR[i*fgNParCh+2]=0.0;
369 fConstraintPR[i*fgNParCh+0]=0.0;
370 fConstraintPR[i*fgNParCh+1]=0.0;
371 fConstraintPR[i*fgNParCh+2]=0.0;
375 void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
376 /// Constrain equation defined by par to value
377 fMillepede->SetGlobalConstraint(par, value);
378 AliInfo("Adding constraint");
381 void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
382 /// Initialize global parameters with par array
383 fMillepede->SetGlobalParameters(par);
384 AliInfo("Init Global Parameters");
387 void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
388 /// Parameter iPar is encourage to vary in [-value;value].
389 /// If value == 0, parameter is fixed
390 fMillepede->SetParSigma(iPar, value);
391 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
394 void AliMUONAlignment::ResetLocalEquation()
396 /// Reset the derivative vectors
397 for(int i=0; i<fNLocal; i++) {
398 fLocalDerivatives[i] = 0.0;
400 for(int i=0; i<fNGlobal; i++) {
401 fGlobalDerivatives[i] = 0.0;
405 void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) {
406 /// Set allowed variation for selected stations based on fDoF and fAllowVar
407 for (Int_t iSt=1; iSt<=5; iSt++) {
408 if (bStOnOff[iSt-1]) {
409 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
410 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
411 for (int i=0; i<fgNParCh; i++) {
412 AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
414 for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
415 FixParameter(j*fgNParCh+i, fAllowVar[i]);
423 void AliMUONAlignment::SetNonLinear(Int_t iPar /* set non linear flag */ ) {
424 /// Set nonlinear flag for parameter iPar
425 fMillepede->SetNonLinear(iPar);
426 AliInfo(Form("Parameter %i set to non linear", iPar));
429 void AliMUONAlignment::LocalEquationX() {
430 /// Define local equation for current track and hit in x coor. measurement
431 // set local derivatives
432 SetLocalDerivative(0, fCosPhi);
433 SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
434 SetLocalDerivative(2, fSinPhi);
435 SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
437 // set global derivatives
438 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -1.);
439 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, 0.);
441 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
442 -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
443 +fCosPhi*(fTrackPos[1]-fDetElemPos[1]));
446 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
447 -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*
448 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
449 +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*
450 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
453 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
456 void AliMUONAlignment::LocalEquationY() {
457 /// Define local equation for current track and hit in y coor. measurement
458 // set local derivatives
459 SetLocalDerivative(0,-fSinPhi);
460 SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
461 SetLocalDerivative(2, fCosPhi);
462 SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
464 // set global derivatives
465 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, 0.);
466 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -1.);
468 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
469 -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
470 -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
473 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
474 -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*
475 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
476 -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*
477 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
480 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
483 void AliMUONAlignment::FillRecPointData() {
484 /// Get information of current hit
485 fClustPos[0] = fRecHit->GetNonBendingCoor();
486 fClustPos[1] = fRecHit->GetBendingCoor();
487 fClustPos[2] = fRecHit->GetZ();
488 fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
489 fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
492 void AliMUONAlignment::FillTrackParamData() {
493 /// Get information of current track at current hit
494 fTrackPos[0] = fTrackParam->GetNonBendingCoor();
495 fTrackPos[1] = fTrackParam->GetBendingCoor();
496 fTrackPos[2] = fTrackParam->GetZ();
497 fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
498 fTrackSlope[1] = fTrackParam->GetBendingSlope();
499 fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
500 fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
503 void AliMUONAlignment::FillDetElemData() {
504 /// Get information of current detection element
505 Double_t lDetElemLocX = 0.;
506 Double_t lDetElemLocY = 0.;
507 Double_t lDetElemLocZ = 0.;
508 fDetElemId = fRecHit->GetDetElemId();
509 fDetElemNumber = fDetElemId%100;
510 for (int iCh=0; iCh<fDetElemId/100-1; iCh++){
511 fDetElemNumber += fgNDetElemCh[iCh];
513 fTransform->Local2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
514 fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
515 if (fDetElemId/100 < 5){
525 void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) {
526 /// Process track; Loop over hits and set local equations
528 // get tclones arrays.
529 fTrackParamAtHit = fTrack->GetTrackParamAtHit();
530 fHitForRecAtHit = fTrack->GetHitForRecAtHit();
532 // get size of arrays
533 Int_t nTrackParam = fTrackParamAtHit->GetEntries();
534 Int_t nHitForRec = fHitForRecAtHit->GetEntries();
535 AliInfo(Form("Number of track param entries : %i ", nTrackParam));
536 AliInfo(Form("Number of hit for rec entries : %i ", nHitForRec));
538 for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
539 fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
540 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
541 if (!fRecHit || !fTrackParam) continue;
542 // fill local variables for this position --> one measurement
545 FillTrackParamData();
546 // if (fDetElemId<500) continue;
547 fTrackPos0[0] = fTrackPos[0];
548 fTrackPos0[1] = fTrackPos[1];
549 fTrackPos0[2] = fTrackPos[2];
550 fTrackSlope0[0] = fTrackSlope[0];
551 fTrackSlope0[1] = fTrackSlope[1];
555 for(Int_t iHit=0; iHit<nHitForRec; iHit++) {
556 // and get new pointers
557 fRecHit = (AliMUONHitForRec *) fHitForRecAtHit->At(iHit);
558 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtHit()->At(iHit);
559 if (!fRecHit || !fTrackParam) continue;
560 // fill local variables for this position --> one measurement
563 FillTrackParamData();
564 // if (fDetElemId<500) continue;
565 AliDebug(1,Form("cluster: %i", iHit));
566 AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
567 AliDebug(1,Form("fDetElemPos[0]: %f\t fDetElemPos[1]: %f\t fDetElemPos[2]: %f\t DetElemID: %i\t ", fDetElemPos[0],fDetElemPos[1],fDetElemPos[2], fDetElemId));
569 AliDebug(1,Form("Track Parameter: %i", iHit));
570 AliDebug(1,Form("x: %f\t y: %f\t z: %f\t slopex: %f\t slopey: %f", fTrackPos[0], fTrackPos[1], fTrackPos[2], fTrackSlope[0], fTrackSlope[1]));
571 AliDebug(1,Form("x0: %f\t y0: %f\t z0: %f\t slopex0: %f\t slopey0: %f", fTrackPos0[0], fTrackPos0[1], fTrackPos0[2], fTrackSlope0[0], fTrackSlope0[1]));
573 fCosPhi = TMath::Cos(fPhi);
574 fSinPhi = TMath::Sin(fPhi);
576 fMeas[0] = fTrackPos[0] - fClustPos[0];
577 fMeas[1] = fTrackPos[1] - fClustPos[1];
580 fMeas[0] = - fClustPos[0];
581 fMeas[1] = - fClustPos[1];
583 AliDebug(1,Form("fMeas[0]: %f\t fMeas[1]: %f\t fSigma[0]: %f\t fSigma[1]: %f", fMeas[0], fMeas[1], fSigma[0], fSigma[1]));
584 // Set local equations
590 void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
591 /// Call local fit for this tracks
592 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
593 if (iRes && !lSingleFit) {
594 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
598 void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
599 /// Call global fit; Global parameters are stored in parameters
600 fMillepede->GlobalFit(parameters,errors,pulls);
601 AliInfo("Done fitting global parameters!");
602 for (int iGlob=0; iGlob<fgNDetElem; iGlob++){
603 printf("%d\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+2]);
607 Double_t AliMUONAlignment::GetParError(Int_t iPar) {
608 /// Get error of parameter iPar
609 Double_t lErr = fMillepede->GetParError(iPar);
613 void AliMUONAlignment::PrintGlobalParameters() {
614 /// Print global parameters
615 fMillepede->PrintGlobalParameters();
618 //_________________________________________________________________________
619 TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *detElemMisAlignment) const
621 /// Realign given transformation by given misalignment and return the misaligned transformation
623 Double_t cartMisAlig[3] = {0,0,0};
624 Double_t angMisAlig[3] = {0,0,0};
625 const Double_t *trans = transform.GetTranslation();
627 // check if the rotation we obtain is not NULL
628 if (transform.GetRotation()) {
629 rot = transform.GetRotation();
632 rot = new TGeoRotation("rot");
633 } // default constructor.
635 cartMisAlig[0] = -detElemMisAlignment[0];
636 cartMisAlig[1] = -detElemMisAlignment[1];
637 angMisAlig[2] = -detElemMisAlignment[2]*180./TMath::Pi();
639 TGeoTranslation newTrans(cartMisAlig[0] + trans[0], cartMisAlig[1] + trans[1], cartMisAlig[2] + trans[2]);
641 rot->RotateX(angMisAlig[0]);
642 rot->RotateY(angMisAlig[1]);
643 rot->RotateZ(angMisAlig[2]);
645 return TGeoCombiTrans(newTrans, *rot);
648 //______________________________________________________________________
649 AliMUONGeometryTransformer *
650 AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
651 double *misAlignments, Bool_t verbose)
654 /////////////////////////////////////////////////////////////////////
655 // Takes the internal geometry module transformers, copies them
656 // and gets the Detection Elements from them.
657 // Takes misalignment parameters and applies these
658 // to the local transform of the Detection Element
659 // Obtains the global transform by multiplying the module transformer
660 // transformation with the local transformation
661 // Applies the global transform to a new detection element
662 // Adds the new detection element to a new module transformer
663 // Adds the new module transformer to a new geometry transformer
664 // Returns the new geometry transformer
666 Double_t lDetElemMisAlignment[3] = {0.,0.,0.};
668 AliMUONGeometryTransformer *newGeometryTransformer =
669 new AliMUONGeometryTransformer(kTRUE);
670 for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++)
671 { // module transformers
673 const AliMUONGeometryModuleTransformer *kModuleTransformer =
674 transformer->GetModuleTransformer(iMt, true);
676 AliMUONGeometryModuleTransformer *newModuleTransformer =
677 new AliMUONGeometryModuleTransformer(iMt);
678 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
680 TGeoCombiTrans moduleTransform =
681 TGeoCombiTrans(*kModuleTransformer->GetTransformation());
682 TGeoCombiTrans *newModuleTransform = new TGeoCombiTrans(moduleTransform);
683 // same module transform as the previous one
684 // no mis align object created
685 newModuleTransformer->SetTransformation(moduleTransform);
687 AliMUONGeometryStore *detElements =
688 kModuleTransformer->GetDetElementStore();
692 ("%i DEs in old GeometryStore %i",
693 detElements->GetNofEntries(), iMt));
694 Int_t iBase = !iMt ? 0 : fgSNDetElemCh[iMt-1];
695 for (Int_t iDe = 0; iDe < detElements->GetNofEntries(); iDe++)
696 { // detection elements.
697 AliMUONGeometryDetElement *detElement =
698 (AliMUONGeometryDetElement *) detElements->GetEntry(iDe);
700 AliFatal("Detection element not found.");
702 /// make a new detection element
703 AliMUONGeometryDetElement *newDetElement =
704 new AliMUONGeometryDetElement(detElement->GetId(),
705 detElement->GetVolumePath());
706 for (int i=0; i<fgNParCh; i++) {
707 AliInfo(Form("iMt %i, iBase %i, iDe %i, iPar %i",iMt, iBase, iDe, (iBase+iDe)*fgNParCh+i));
708 lDetElemMisAlignment[i] =
709 (iMt<fgNCh) ? misAlignments[(iBase+iDe)*fgNParCh+i] : 0.;
711 // local transformation of this detection element.
712 TGeoCombiTrans localTransform
713 = TGeoCombiTrans(*detElement->GetLocalTransformation());
714 TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
715 newDetElement->SetLocalTransformation(newLocalTransform);
717 // global transformation
718 TGeoHMatrix newGlobalTransform =
719 AliMUONGeometryBuilder::Multiply(*newModuleTransform,
721 newDetElement->SetGlobalTransformation(newGlobalTransform);
723 // add this det element to module
724 newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
726 // Get delta transformation:
727 // Tdelta = Tnew * Told.inverse
728 TGeoHMatrix deltaTransform
729 = AliMUONGeometryBuilder::Multiply(
731 detElement->GetGlobalTransformation()->Inverse());
733 // Create mis alignment matrix
734 newGeometryTransformer
735 ->AddMisAlignDetElement(detElement->GetId(), deltaTransform);
738 AliInfo(Form("Added module transformer %i to the transformer", iMt));
739 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
741 return newGeometryTransformer;