]> 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
AliEveTrackCounter, AliEveTrackCounterEditor
[u/mrichter/AliRoot.git] / MUON / AliMUONAlignment.cxx
CommitLineData
010eb601 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//-----------------------------------------------------------------------------
c25af45b 19/// \class AliMUONAlignment
20/// Alignment class for the ALICE DiMuon spectrometer
010eb601 21///
22/// MUON specific alignment class which interface to AliMillepede.
23/// For each track ProcessTrack calculates the local and global derivatives
96ebe67e 24/// at each cluster and fill the corresponding local equations. Provide methods
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"
010eb601 38#include "AliMUONConstants.h"
010eb601 39#include "AliMillepede.h"
40
8c95578a 41#include "AliMpExMap.h"
630711ed 42#include "AliMpExMapIterator.h"
8c95578a 43
4818a9b7 44#include "AliAlignObjMatrix.h"
8c95578a 45#include "AliLog.h"
46
2060b217 47#include "TMath.h"
4818a9b7 48#include "TMatrixDSym.h"
8c95578a 49#include "TSystem.h"
50
78649106 51/// \cond CLASSIMP
010eb601 52ClassImp(AliMUONAlignment)
78649106 53/// \endcond
54
010eb601 55 Int_t AliMUONAlignment::fgNDetElem = 4*2+4*2+18*2+26*2+26*2;
56 Int_t AliMUONAlignment::fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
57 Int_t AliMUONAlignment::fgSNDetElemCh[10] = {4,8,12,16,34,52,78,104,130,156};
58 Int_t AliMUONAlignment::fgNParCh = 3;
3b2890be 59 Int_t AliMUONAlignment::fgNTrkMod = 16;
010eb601 60 Int_t AliMUONAlignment::fgNCh = 10;
61 Int_t AliMUONAlignment::fgNSt = 5;
62
8395bff1 63AliMUONAlignment::AliMUONAlignment()
64 : TObject(),
65 fBFieldOn(kTRUE),
66 fStartFac(16.),
67 fResCutInitial(100.),
68 fResCut(100.),
69 fMillepede(0),
96ebe67e 70 fTrackParamAtCluster(0),
8395bff1 71 fTrack(0),
96ebe67e 72 fCluster(0),
8395bff1 73 fTrackParam(0),
74 fNGlobal(fgNDetElem*fgNParCh),
75 fNLocal(4),
76 fNStdDev(3),
77 fDetElemId(0),
78 fDetElemNumber(0),
79 fPhi(0.),
80 fCosPhi(1.),
81 fSinPhi(0.),
82 fTransform(0)
010eb601 83{
84 /// Default constructor, setting define alignment parameters
85 fSigma[0] = 1.0e-1;
86 fSigma[1] = 1.0e-2;
010eb601 87
88 fDoF[0] = kTRUE; fDoF[1] = kTRUE; fDoF[2] = kTRUE;
89 fAllowVar[0] = 0.05; fAllowVar[1] = 0.05; fAllowVar[2] = 0.001;
8395bff1 90
010eb601 91 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));
92
93 fMillepede = new AliMillepede();
94
95 Init(fNGlobal, fNLocal, fNStdDev);
96
97 ResetLocalEquation();
98 AliInfo("Parameters initialized to zero");
99
100}
101
102AliMUONAlignment::~AliMUONAlignment() {
103 /// Destructor
104}
105
106void AliMUONAlignment::Init(Int_t nGlobal, /* number of global paramers */
107 Int_t nLocal, /* number of local parameters */
108 Int_t nStdDev /* std dev cut */ )
109{
110 /// Initialization of AliMillepede. Fix parameters, define constraints ...
111 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
112
3b2890be 113 Bool_t bStOnOff[5] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
114 Bool_t bSpecLROnOff[2] = {kTRUE,kTRUE};
010eb601 115
116 AllowVariations(bStOnOff);
117
118 // Fix parameters or add constraints here
119 for (Int_t iSt=0; iSt<5; iSt++)
120 if (!bStOnOff[iSt]) FixStation(iSt+1);
121
3b2890be 122 FixHalfSpectrometer(bStOnOff,bSpecLROnOff);
123
010eb601 124 ResetConstraints();
3b2890be 125
126 // Define global constrains to be applied
127 // X, Y, P, XvsZ, YvsZ, PvsZ, XvsY, YvsY, PvsY
128 Bool_t bVarXYT[9] = {kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
129 Bool_t bDetTLBR[4] = {kFALSE,kTRUE,kFALSE,kTRUE};
130 // AddConstraints(bStOnOff,bVarXYT,bDetTLBR,bSpecLROnOff);
131
132 // Other possible way to add constrains
010eb601 133 bVarXYT[0] = kFALSE; bVarXYT[1] = kFALSE; bVarXYT[2] = kTRUE;
134 bDetTLBR[0] = kFALSE; bDetTLBR[1] = kTRUE; bDetTLBR[2] = kFALSE; bDetTLBR[3] = kFALSE;
3b2890be 135// AddConstraints(bStOnOff,bVarXYT,bDetTLBR);
136
137 bVarXYT[0] = kTRUE; bVarXYT[1] = kTRUE; bVarXYT[2] = kFALSE;
138 // AddConstraints(bStOnOff,bVarXYT);
010eb601 139
140 // Set iterations
141 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
142}
143
144void AliMUONAlignment::FixStation(Int_t iSt){
145 /// Fix all detection elements of station iSt
146 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
147 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
148 for (Int_t i = iDetElemFirst; i < iDetElemLast; i++){
149 FixParameter(i*fgNParCh+0, 0.0);
150 FixParameter(i*fgNParCh+1, 0.0);
151 FixParameter(i*fgNParCh+2, 0.0);
152 }
153}
154
3b2890be 155void AliMUONAlignment::FixHalfSpectrometer(Bool_t *lStOnOff, Bool_t *lSpecLROnOff){
156 /// Fix left or right detector
157 for (Int_t i = 0; i < fgNDetElem; i++){
158 Int_t iCh=0;
159 for (iCh=1; iCh<=fgNCh; iCh++){
160 if (i<fgSNDetElemCh[iCh-1]) break;
161 }
162 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
163 if (iSt){
164 Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
165 if (iCh>=1 && iCh<=4){
166 if ((lDetElemNumber==1 || lDetElemNumber==2) && !lSpecLROnOff[0]){ // From track crossings
167 FixParameter(i*fgNParCh+0, 0.0);
168 FixParameter(i*fgNParCh+1, 0.0);
169 FixParameter(i*fgNParCh+2, 0.0);
170 }
171 if ((lDetElemNumber==0 || lDetElemNumber==3) && !lSpecLROnOff[1]){ // From track crossings
172 FixParameter(i*fgNParCh+0, 0.0);
173 FixParameter(i*fgNParCh+1, 0.0);
174 FixParameter(i*fgNParCh+2, 0.0);
175 }
176 }
177 if (iCh>=5 && iCh<=6){
178 if ((lDetElemNumber>=5&&lDetElemNumber<=13) && !lSpecLROnOff[0]){
179 FixParameter(i*fgNParCh+0, 0.0);
180 FixParameter(i*fgNParCh+1, 0.0);
181 FixParameter(i*fgNParCh+2, 0.0);
182 }
183 if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
184 (lDetElemNumber>=14&&lDetElemNumber<=17)) && !lSpecLROnOff[1]){
185 FixParameter(i*fgNParCh+0, 0.0);
186 FixParameter(i*fgNParCh+1, 0.0);
187 FixParameter(i*fgNParCh+2, 0.0);
188 }
189 }
190 if (iCh>=7 && iCh<=10){
191 if ((lDetElemNumber>=7&&lDetElemNumber<=19) && !lSpecLROnOff[0]){
192 FixParameter(i*fgNParCh+0, 0.0);
193 FixParameter(i*fgNParCh+1, 0.0);
194 FixParameter(i*fgNParCh+2, 0.0);
195 }
196 if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
197 (lDetElemNumber>=20&&lDetElemNumber<=25)) && !lSpecLROnOff[1]){
198 FixParameter(i*fgNParCh+0, 0.0);
199 FixParameter(i*fgNParCh+1, 0.0);
200 FixParameter(i*fgNParCh+2, 0.0);
201 }
202 }
203 }
204 }
205}
206
010eb601 207void AliMUONAlignment::SetNonLinear(Bool_t *lStOnOff,Bool_t *lVarXYT){
208 /// Set non linear parameter flag selected stations and degrees of freedom
209 for (Int_t i = 0; i < fgNDetElem; i++){
210 Int_t iCh=0;
211 for (iCh=1; iCh<=fgNCh; iCh++){
212 if (i<fgSNDetElemCh[iCh-1]) break;
213 }
214 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
215 if (iSt){
216 if (lVarXYT[0]) { // X constraint
217 SetNonLinear(i*fgNParCh+0);
218 }
219 if (lVarXYT[1]) { // Y constraint
220 SetNonLinear(i*fgNParCh+1);
221 }
222 if (lVarXYT[2]) { // T constraint
223 SetNonLinear(i*fgNParCh+2);
224 }
225 }
226 }
227}
228
229void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT){
230 /// Add constraint equations for selected stations and degrees of freedom
231 for (Int_t i = 0; i < fgNDetElem; i++){
232 Int_t iCh=0;
233 for (iCh=1; iCh<=fgNCh; iCh++){
234 if (i<fgSNDetElemCh[iCh-1]) break;
235 }
236 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
237 if (iSt){
238 if (lVarXYT[0]) { // X constraint
239 fConstraintX[i*fgNParCh+0]=1.0;
240 }
241 if (lVarXYT[1]) { // Y constraint
242 fConstraintY[i*fgNParCh+1]=1.0;
243 }
244 if (lVarXYT[2]) { // T constraint
245 fConstraintP[i*fgNParCh+2]=1.0;
246 }
247 }
248 }
249 if (lVarXYT[0]) { // X constraint
250 AddConstraint(fConstraintX,0.0);
251 }
252 if (lVarXYT[1]) { // Y constraint
253 AddConstraint(fConstraintY,0.0);
254 }
255 if (lVarXYT[2]) { // T constraint
256 AddConstraint(fConstraintP,0.0);
257 }
258}
259
3b2890be 260void AliMUONAlignment::AddConstraints(Bool_t *lStOnOff,Bool_t *lVarXYT, Bool_t *lDetTLBR, Bool_t *lSpecLROnOff){
261 /// Add constraint equations for selected stations, degrees of freedom and detector half
262 Double_t lDetElemLocX = 0.;
263 Double_t lDetElemLocY = 0.;
264 Double_t lDetElemLocZ = 0.;
265 Double_t lDetElemGloX = 0.;
266 Double_t lDetElemGloY = 0.;
267 Double_t lDetElemGloZ = 0.;
268 Double_t lMeanY = 0.;
269 Double_t lSigmaY = 0.;
270 Double_t lMeanZ = 0.;
271 Double_t lSigmaZ = 0.;
272 Int_t lNDetElem = 0;
010eb601 273 for (Int_t i = 0; i < fgNDetElem; i++){
274 Int_t iCh=0;
275 for (iCh=1; iCh<=fgNCh; iCh++){
276 if (i<fgSNDetElemCh[iCh-1]) break;
277 }
278 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
279 if (iSt){
3b2890be 280 Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
281 Int_t lDetElemId = iCh*100+lDetElemNumber;
282 fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
283 lDetElemGloX,lDetElemGloY,lDetElemGloZ);
284 if (iCh>=1 && iCh<=4){
285 if ((lDetElemNumber==1 || lDetElemNumber==2) && lSpecLROnOff[0]){ // From track crossings
286 lMeanY += lDetElemGloY;
287 lSigmaY += lDetElemGloY*lDetElemGloY;
288 lMeanZ += lDetElemGloZ;
289 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
290 lNDetElem++;
291 }
292 if ((lDetElemNumber==0 || lDetElemNumber==3) && lSpecLROnOff[1]){ // From track crossings
293 lMeanY += lDetElemGloY;
294 lSigmaY += lDetElemGloY*lDetElemGloY;
295 lMeanZ += lDetElemGloZ;
296 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
297 lNDetElem++;
298 }
299 }
300 if (iCh>=5 && iCh<=6){
301 if ((lDetElemNumber>=5&&lDetElemNumber<=13) && lSpecLROnOff[0]){
302 lMeanY += lDetElemGloY;
303 lSigmaY += lDetElemGloY*lDetElemGloY;
304 lMeanZ += lDetElemGloZ;
305 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
306 lNDetElem++;
307 }
308 if (((lDetElemNumber>=0&&lDetElemNumber<=4) ||
309 (lDetElemNumber>=14&&lDetElemNumber<=17)) && lSpecLROnOff[1]){
310 lMeanY += lDetElemGloY;
311 lSigmaY += lDetElemGloY*lDetElemGloY;
312 lMeanZ += lDetElemGloZ;
313 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
314 lNDetElem++;
315 }
316 }
317 if (iCh>=7 && iCh<=10){
318 if ((lDetElemNumber>=7&&lDetElemNumber<=19) && lSpecLROnOff[0]){
319 lMeanY += lDetElemGloY;
320 lSigmaY += lDetElemGloY*lDetElemGloY;
321 lMeanZ += lDetElemGloZ;
322 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
323 lNDetElem++;
324 }
325 if (((lDetElemNumber>=0&&lDetElemNumber<=6) ||
326 (lDetElemNumber>=20&&lDetElemNumber<=25)) && lSpecLROnOff[1]){
327 lMeanY += lDetElemGloY;
328 lSigmaY += lDetElemGloY*lDetElemGloY;
329 lMeanZ += lDetElemGloZ;
330 lSigmaZ += lDetElemGloZ*lDetElemGloZ;
331 lNDetElem++;
332 }
333 }
334 }
335 }
336 lMeanY /= lNDetElem;
337 lSigmaY /= lNDetElem;
338 lSigmaY = TMath::Sqrt(lSigmaY-lMeanY*lMeanY);
339 lMeanZ /= lNDetElem;
340 lSigmaZ /= lNDetElem;
341 lSigmaZ = TMath::Sqrt(lSigmaZ-lMeanZ*lMeanZ);
342 AliInfo(Form("Used %i DetElem, MeanZ= %f , SigmaZ= %f", lNDetElem,lMeanZ,lSigmaZ));
343
344 for (Int_t i = 0; i < fgNDetElem; i++){
345 Int_t iCh=0;
346 for (iCh=1; iCh<=fgNCh; iCh++){
347 if (i<fgSNDetElemCh[iCh-1]) break;
348 }
349 Int_t iSt = lStOnOff[(iCh-1)/2] ? (iCh+1)/2 : 0;
350 if (iSt){
351 Int_t lDetElemNumber = (iCh==1) ? i : i-fgSNDetElemCh[iCh-2];
352 Int_t lDetElemId = iCh*100+lDetElemNumber;
353 fTransform->Local2Global(lDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
354 lDetElemGloX,lDetElemGloY,lDetElemGloZ);
010eb601 355 if (lVarXYT[0]) { // X constraint
356 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXT,0); // Top half
357 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXL,0); // Left half
358 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXB,0); // Bottom half
359 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXR,0); // Right half
360 }
3b2890be 361 if (lVarXYT[1]) { // Y constraint
010eb601 362 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYT,1); // Top half
363 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYL,1); // Left half
364 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYB,1); // Bottom half
365 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYR,1); // Right half
366 }
3b2890be 367 if (lVarXYT[2]) { // P constraint
010eb601 368 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPT,2); // Top half
369 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPL,2); // Left half
370 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPB,2); // Bottom half
371 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPR,2); // Right half
372 }
3b2890be 373 if (lVarXYT[3]) { // X-Z shearing
374 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXZT,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
375 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXZL,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
376 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXZB,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
377 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXZR,0,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
378 }
379 if (lVarXYT[4]) { // Y-Z shearing
380 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYZT,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
381 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYZL,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
382 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYZB,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
383 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYZR,1,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
384 }
385 if (lVarXYT[5]) { // P-Z rotation
386 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPZT,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Top half
387 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPZL,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Left half
388 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPZB,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Bottom half
389 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPZR,2,(lDetElemGloZ-lMeanZ)/lSigmaZ); // Right half
390 }
391 if (lVarXYT[6]) { // X-Y shearing
392 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintXYT,0,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
393 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintXYL,0,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
394 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintXYB,0,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
395 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintXYR,0,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
396 }
397 if (lVarXYT[7]) { // Y-Y scaling
398 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintYYT,1,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
399 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintYYL,1,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
400 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintYYB,1,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
401 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintYYR,1,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
402 }
403 if (lVarXYT[8]) { // P-Y rotation
404 if (lDetTLBR[0]) ConstrainT(i,iCh,fConstraintPYT,2,(lDetElemGloY-lMeanY)/lSigmaY); // Top half
405 if (lDetTLBR[1]) ConstrainL(i,iCh,fConstraintPYL,2,(lDetElemGloY-lMeanY)/lSigmaY); // Left half
406 if (lDetTLBR[2]) ConstrainB(i,iCh,fConstraintPYB,2,(lDetElemGloY-lMeanY)/lSigmaY); // Bottom half
407 if (lDetTLBR[3]) ConstrainR(i,iCh,fConstraintPYR,2,(lDetElemGloY-lMeanY)/lSigmaY); // Right half
408 }
010eb601 409 }
410 }
411 if (lVarXYT[0]) { // X constraint
412 if (lDetTLBR[0]) AddConstraint(fConstraintXT,0.0); // Top half
413 if (lDetTLBR[1]) AddConstraint(fConstraintXL,0.0); // Left half
414 if (lDetTLBR[2]) AddConstraint(fConstraintXB,0.0); // Bottom half
415 if (lDetTLBR[3]) AddConstraint(fConstraintXR,0.0); // Right half
416 }
417 if (lVarXYT[1]) { // Y constraint
418 if (lDetTLBR[0]) AddConstraint(fConstraintYT,0.0); // Top half
419 if (lDetTLBR[1]) AddConstraint(fConstraintYL,0.0); // Left half
420 if (lDetTLBR[2]) AddConstraint(fConstraintYB,0.0); // Bottom half
421 if (lDetTLBR[3]) AddConstraint(fConstraintYR,0.0); // Right half
422 }
423 if (lVarXYT[2]) { // T constraint
424 if (lDetTLBR[0]) AddConstraint(fConstraintPT,0.0); // Top half
425 if (lDetTLBR[1]) AddConstraint(fConstraintPL,0.0); // Left half
426 if (lDetTLBR[2]) AddConstraint(fConstraintPB,0.0); // Bottom half
427 if (lDetTLBR[3]) AddConstraint(fConstraintPR,0.0); // Right half
428 }
3b2890be 429 if (lVarXYT[3]) { // X-Z constraint
430 if (lDetTLBR[0]) AddConstraint(fConstraintXZT,0.0); // Top half
431 if (lDetTLBR[1]) AddConstraint(fConstraintXZL,0.0); // Left half
432 if (lDetTLBR[2]) AddConstraint(fConstraintXZB,0.0); // Bottom half
433 if (lDetTLBR[3]) AddConstraint(fConstraintXZR,0.0); // Right half
434 }
435 if (lVarXYT[4]) { // Y-Z constraint
436 if (lDetTLBR[0]) AddConstraint(fConstraintYZT,0.0); // Top half
437 if (lDetTLBR[1]) AddConstraint(fConstraintYZL,0.0); // Left half
438 if (lDetTLBR[2]) AddConstraint(fConstraintYZB,0.0); // Bottom half
439 if (lDetTLBR[3]) AddConstraint(fConstraintYZR,0.0); // Right half
440 }
441 if (lVarXYT[5]) { // P-Z constraint
442 if (lDetTLBR[0]) AddConstraint(fConstraintPZT,0.0); // Top half
443 if (lDetTLBR[1]) AddConstraint(fConstraintPZL,0.0); // Left half
444 if (lDetTLBR[2]) AddConstraint(fConstraintPZB,0.0); // Bottom half
445 if (lDetTLBR[3]) AddConstraint(fConstraintPZR,0.0); // Right half
446 }
447 if (lVarXYT[6]) { // X-Y constraint
448 if (lDetTLBR[0]) AddConstraint(fConstraintXYT,0.0); // Top half
449 if (lDetTLBR[1]) AddConstraint(fConstraintXYL,0.0); // Left half
450 if (lDetTLBR[2]) AddConstraint(fConstraintXYB,0.0); // Bottom half
451 if (lDetTLBR[3]) AddConstraint(fConstraintXYR,0.0); // Right half
452 }
453 if (lVarXYT[7]) { // Y-Y constraint
454 if (lDetTLBR[0]) AddConstraint(fConstraintYYT,0.0); // Top half
455 if (lDetTLBR[1]) AddConstraint(fConstraintYYL,0.0); // Left half
456 if (lDetTLBR[2]) AddConstraint(fConstraintYYB,0.0); // Bottom half
457 if (lDetTLBR[3]) AddConstraint(fConstraintYYR,0.0); // Right half
458 }
459 if (lVarXYT[8]) { // P-Y constraint
460 if (lDetTLBR[0]) AddConstraint(fConstraintPYT,0.0); // Top half
461 if (lDetTLBR[1]) AddConstraint(fConstraintPYL,0.0); // Left half
462 if (lDetTLBR[2]) AddConstraint(fConstraintPYB,0.0); // Bottom half
463 if (lDetTLBR[3]) AddConstraint(fConstraintPYR,0.0); // Right half
464 }
010eb601 465}
466
8ec2f43d 467void AliMUONAlignment::ConstrainT(Int_t lDetElem, Int_t lCh, Double_t *lConstraintT, Int_t iVar, Double_t /*lWeight*/){
010eb601 468 /// Set constrain equation for top half of spectrometer
469 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
470 if (lCh>=1 && lCh<=4){
471 if (lDetElemNumber==0 || lDetElemNumber==1){ // From track crossings
472 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
473 }
474 }
475 if (lCh>=5 && lCh<=6){
476 if (lDetElemNumber>=0&&lDetElemNumber<=9){
477 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
478 }
479 }
480 if (lCh>=7 && lCh<=10){
481 if (lDetElemNumber>=0&&lDetElemNumber<=13){
482 lConstraintT[lDetElem*fgNParCh+iVar]=1.0;
483 }
484 }
485}
486
3b2890be 487void AliMUONAlignment::ConstrainL(Int_t lDetElem, Int_t lCh, Double_t *lConstraintL, Int_t iVar, Double_t lWeight){
010eb601 488 /// Set constrain equation for left half of spectrometer
489 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
490 if (lCh>=1 && lCh<=4){
491 if (lDetElemNumber==1 || lDetElemNumber==2){ // From track crossings
3b2890be 492 lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 493 }
494 }
495 if (lCh>=5 && lCh<=6){
496 if (lDetElemNumber>=5&&lDetElemNumber<=13){
3b2890be 497 lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 498 }
499 }
500 if (lCh>=7 && lCh<=10){
501 if (lDetElemNumber>=7&&lDetElemNumber<=19){
3b2890be 502 lConstraintL[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 503 }
504 }
505}
506
8ec2f43d 507void AliMUONAlignment::ConstrainB(Int_t lDetElem, Int_t lCh, Double_t *lConstraintB, Int_t iVar, Double_t /*lWeight*/){
010eb601 508 /// Set constrain equation for bottom half of spectrometer
509 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
510 if (lCh>=1 && lCh<=4){
3b2890be 511 if (lDetElemNumber==2 || lDetElemNumber==3){ // From track crossings
010eb601 512 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
513 }
514 }
515 if (lCh>=5 && lCh<=6){
516 if ((lDetElemNumber>=9&&lDetElemNumber<=17) ||
517 (lDetElemNumber==0)){
518 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
519 }
520 }
521 if (lCh>=7 && lCh<=10){
522 if ((lDetElemNumber>=13&&lDetElemNumber<=25) ||
523 (lDetElemNumber==0)){
524 lConstraintB[lDetElem*fgNParCh+iVar]=1.0;
525 }
526 }
527}
528
3b2890be 529void AliMUONAlignment::ConstrainR(Int_t lDetElem, Int_t lCh, Double_t *lConstraintR, Int_t iVar, Double_t lWeight){
010eb601 530 /// Set constrain equation for right half of spectrometer
531 Int_t lDetElemNumber = (lCh==1) ? lDetElem : lDetElem-fgSNDetElemCh[lCh-2];
532 if (lCh>=1 && lCh<=4){
3b2890be 533 if (lDetElemNumber==0 || lDetElemNumber==3){ // From track crossings
534 lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 535 }
536 }
537 if (lCh>=5 && lCh<=6){
538 if ((lDetElemNumber>=0&&lDetElemNumber<=4) ||
539 (lDetElemNumber>=14&&lDetElemNumber<=17)){
3b2890be 540 lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 541 }
542 }
543 if (lCh>=7 && lCh<=10){
544 if ((lDetElemNumber>=0&&lDetElemNumber<=6) ||
545 (lDetElemNumber>=20&&lDetElemNumber<=25)){
3b2890be 546 lConstraintR[lDetElem*fgNParCh+iVar]=lWeight;
010eb601 547 }
548 }
549}
550
551void AliMUONAlignment::ResetConstraints(){
552 /// Reset all constraint equations
553 for (Int_t i = 0; i < fgNDetElem; i++){
3b2890be 554 fConstraintX[i*fgNParCh+0]=0.0;
555 fConstraintX[i*fgNParCh+1]=0.0;
556 fConstraintX[i*fgNParCh+2]=0.0;
557 fConstraintY[i*fgNParCh+0]=0.0;
558 fConstraintY[i*fgNParCh+1]=0.0;
559 fConstraintY[i*fgNParCh+2]=0.0;
560 fConstraintP[i*fgNParCh+0]=0.0;
561 fConstraintP[i*fgNParCh+1]=0.0;
562 fConstraintP[i*fgNParCh+2]=0.0;
563 fConstraintXT[i*fgNParCh+0]=0.0;
564 fConstraintXT[i*fgNParCh+1]=0.0;
565 fConstraintXT[i*fgNParCh+2]=0.0;
566 fConstraintYT[i*fgNParCh+0]=0.0;
567 fConstraintYT[i*fgNParCh+1]=0.0;
568 fConstraintYT[i*fgNParCh+2]=0.0;
569 fConstraintPT[i*fgNParCh+0]=0.0;
570 fConstraintPT[i*fgNParCh+1]=0.0;
571 fConstraintPT[i*fgNParCh+2]=0.0;
572 fConstraintXZT[i*fgNParCh+0]=0.0;
573 fConstraintXZT[i*fgNParCh+1]=0.0;
574 fConstraintXZT[i*fgNParCh+2]=0.0;
575 fConstraintYZT[i*fgNParCh+0]=0.0;
576 fConstraintYZT[i*fgNParCh+1]=0.0;
577 fConstraintYZT[i*fgNParCh+2]=0.0;
578 fConstraintPZT[i*fgNParCh+0]=0.0;
579 fConstraintPZT[i*fgNParCh+1]=0.0;
580 fConstraintPZT[i*fgNParCh+2]=0.0;
581 fConstraintXYT[i*fgNParCh+0]=0.0;
582 fConstraintXYT[i*fgNParCh+1]=0.0;
583 fConstraintXYT[i*fgNParCh+2]=0.0;
584 fConstraintYYT[i*fgNParCh+0]=0.0;
585 fConstraintYYT[i*fgNParCh+1]=0.0;
586 fConstraintYYT[i*fgNParCh+2]=0.0;
587 fConstraintPYT[i*fgNParCh+0]=0.0;
588 fConstraintPYT[i*fgNParCh+1]=0.0;
589 fConstraintPYT[i*fgNParCh+2]=0.0;
590 fConstraintXL[i*fgNParCh+0]=0.0;
591 fConstraintXL[i*fgNParCh+1]=0.0;
592 fConstraintXL[i*fgNParCh+2]=0.0;
593 fConstraintYL[i*fgNParCh+0]=0.0;
594 fConstraintYL[i*fgNParCh+1]=0.0;
595 fConstraintYL[i*fgNParCh+2]=0.0;
596 fConstraintPL[i*fgNParCh+0]=0.0;
597 fConstraintPL[i*fgNParCh+1]=0.0;
598 fConstraintPL[i*fgNParCh+2]=0.0;
599 fConstraintXZL[i*fgNParCh+0]=0.0;
600 fConstraintXZL[i*fgNParCh+1]=0.0;
601 fConstraintXZL[i*fgNParCh+2]=0.0;
602 fConstraintYZL[i*fgNParCh+0]=0.0;
603 fConstraintYZL[i*fgNParCh+1]=0.0;
604 fConstraintYZL[i*fgNParCh+2]=0.0;
605 fConstraintPZL[i*fgNParCh+0]=0.0;
606 fConstraintPZL[i*fgNParCh+1]=0.0;
607 fConstraintPZL[i*fgNParCh+2]=0.0;
608 fConstraintXYL[i*fgNParCh+0]=0.0;
609 fConstraintXYL[i*fgNParCh+1]=0.0;
610 fConstraintXYL[i*fgNParCh+2]=0.0;
611 fConstraintYYL[i*fgNParCh+0]=0.0;
612 fConstraintYYL[i*fgNParCh+1]=0.0;
613 fConstraintYYL[i*fgNParCh+2]=0.0;
614 fConstraintPYL[i*fgNParCh+0]=0.0;
615 fConstraintPYL[i*fgNParCh+1]=0.0;
616 fConstraintPYL[i*fgNParCh+2]=0.0;
617 fConstraintXB[i*fgNParCh+0]=0.0;
618 fConstraintXB[i*fgNParCh+1]=0.0;
619 fConstraintXB[i*fgNParCh+2]=0.0;
620 fConstraintYB[i*fgNParCh+0]=0.0;
621 fConstraintYB[i*fgNParCh+1]=0.0;
622 fConstraintYB[i*fgNParCh+2]=0.0;
623 fConstraintPB[i*fgNParCh+0]=0.0;
624 fConstraintPB[i*fgNParCh+1]=0.0;
625 fConstraintPB[i*fgNParCh+2]=0.0;
626 fConstraintXZB[i*fgNParCh+0]=0.0;
627 fConstraintXZB[i*fgNParCh+1]=0.0;
628 fConstraintXZB[i*fgNParCh+2]=0.0;
629 fConstraintYZB[i*fgNParCh+0]=0.0;
630 fConstraintYZB[i*fgNParCh+1]=0.0;
631 fConstraintYZB[i*fgNParCh+2]=0.0;
632 fConstraintPZB[i*fgNParCh+0]=0.0;
633 fConstraintPZB[i*fgNParCh+1]=0.0;
634 fConstraintPZB[i*fgNParCh+2]=0.0;
635 fConstraintXYB[i*fgNParCh+0]=0.0;
636 fConstraintXYB[i*fgNParCh+1]=0.0;
637 fConstraintXYB[i*fgNParCh+2]=0.0;
638 fConstraintYYB[i*fgNParCh+0]=0.0;
639 fConstraintYYB[i*fgNParCh+1]=0.0;
640 fConstraintYYB[i*fgNParCh+2]=0.0;
641 fConstraintPYB[i*fgNParCh+0]=0.0;
642 fConstraintPYB[i*fgNParCh+1]=0.0;
643 fConstraintPYB[i*fgNParCh+2]=0.0;
644 fConstraintXR[i*fgNParCh+0]=0.0;
645 fConstraintXR[i*fgNParCh+1]=0.0;
646 fConstraintXR[i*fgNParCh+2]=0.0;
647 fConstraintYR[i*fgNParCh+0]=0.0;
648 fConstraintYR[i*fgNParCh+1]=0.0;
649 fConstraintYR[i*fgNParCh+2]=0.0;
650 fConstraintPR[i*fgNParCh+0]=0.0;
651 fConstraintPR[i*fgNParCh+1]=0.0;
652 fConstraintPR[i*fgNParCh+2]=0.0;
653 fConstraintXZR[i*fgNParCh+0]=0.0;
654 fConstraintXZR[i*fgNParCh+1]=0.0;
655 fConstraintXZR[i*fgNParCh+2]=0.0;
656 fConstraintYZR[i*fgNParCh+0]=0.0;
657 fConstraintYZR[i*fgNParCh+1]=0.0;
658 fConstraintYZR[i*fgNParCh+2]=0.0;
659 fConstraintPZR[i*fgNParCh+0]=0.0;
660 fConstraintPZR[i*fgNParCh+1]=0.0;
661 fConstraintPZR[i*fgNParCh+2]=0.0;
662 fConstraintPZR[i*fgNParCh+0]=0.0;
663 fConstraintPZR[i*fgNParCh+1]=0.0;
664 fConstraintPZR[i*fgNParCh+2]=0.0;
665 fConstraintXYR[i*fgNParCh+0]=0.0;
666 fConstraintXYR[i*fgNParCh+1]=0.0;
667 fConstraintXYR[i*fgNParCh+2]=0.0;
668 fConstraintYYR[i*fgNParCh+0]=0.0;
669 fConstraintYYR[i*fgNParCh+1]=0.0;
670 fConstraintYYR[i*fgNParCh+2]=0.0;
671 fConstraintPYR[i*fgNParCh+0]=0.0;
672 fConstraintPYR[i*fgNParCh+1]=0.0;
673 fConstraintPYR[i*fgNParCh+2]=0.0;
010eb601 674 }
675}
676
677void AliMUONAlignment::AddConstraint(Double_t *par, Double_t value) {
678 /// Constrain equation defined by par to value
679 fMillepede->SetGlobalConstraint(par, value);
680 AliInfo("Adding constraint");
681}
682
683void AliMUONAlignment::InitGlobalParameters(Double_t *par) {
684 /// Initialize global parameters with par array
685 fMillepede->SetGlobalParameters(par);
686 AliInfo("Init Global Parameters");
687}
688
689void AliMUONAlignment::FixParameter(Int_t iPar, Double_t value) {
690 /// Parameter iPar is encourage to vary in [-value;value].
691 /// If value == 0, parameter is fixed
692 fMillepede->SetParSigma(iPar, value);
693 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
694}
695
696void AliMUONAlignment::ResetLocalEquation()
697{
698 /// Reset the derivative vectors
699 for(int i=0; i<fNLocal; i++) {
700 fLocalDerivatives[i] = 0.0;
701 }
702 for(int i=0; i<fNGlobal; i++) {
703 fGlobalDerivatives[i] = 0.0;
704 }
705}
706
707void AliMUONAlignment::AllowVariations(Bool_t *bStOnOff) {
708 /// Set allowed variation for selected stations based on fDoF and fAllowVar
709 for (Int_t iSt=1; iSt<=5; iSt++) {
710 if (bStOnOff[iSt-1]) {
711 Int_t iDetElemFirst = (iSt>1) ? fgSNDetElemCh[2*(iSt-1)-1] : 0;
712 Int_t iDetElemLast = fgSNDetElemCh[2*(iSt)-1];
713 for (int i=0; i<fgNParCh; i++) {
714 AliDebug(1,Form("fDoF[%d]= %d",i,fDoF[i]));
715 if (fDoF[i]) {
716 for (Int_t j=iDetElemFirst; j<iDetElemLast; j++){
717 FixParameter(j*fgNParCh+i, fAllowVar[i]);
718 }
719 }
720 }
721 }
722 }
723}
724
725void AliMUONAlignment::SetNonLinear(Int_t iPar /* set non linear flag */ ) {
726 /// Set nonlinear flag for parameter iPar
727 fMillepede->SetNonLinear(iPar);
728 AliInfo(Form("Parameter %i set to non linear", iPar));
729}
730
731void AliMUONAlignment::LocalEquationX() {
96ebe67e 732 /// Define local equation for current track and cluster in x coor. measurement
010eb601 733 // set local derivatives
734 SetLocalDerivative(0, fCosPhi);
735 SetLocalDerivative(1, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
736 SetLocalDerivative(2, fSinPhi);
737 SetLocalDerivative(3, fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
738
739 // set global derivatives
3b2890be 740 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, -fCosPhi);
741 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fSinPhi);
010eb601 742 if (fBFieldOn){
743 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
744 -fSinPhi*(fTrackPos[0]-fDetElemPos[0])
745 +fCosPhi*(fTrackPos[1]-fDetElemPos[1]));
746 }
747 else {
748 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
749 -fSinPhi*(fTrackPos0[0]+fTrackSlope0[0]*
750 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
751 +fCosPhi*(fTrackPos0[1]+fTrackSlope0[1]*
752 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
753 }
754
755 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[0], fSigma[0]);
756}
757
758void AliMUONAlignment::LocalEquationY() {
96ebe67e 759 /// Define local equation for current track and cluster in y coor. measurement
010eb601 760 // set local derivatives
761 SetLocalDerivative(0,-fSinPhi);
762 SetLocalDerivative(1,-fSinPhi * (fTrackPos[2] - fTrackPos0[2]));
763 SetLocalDerivative(2, fCosPhi);
764 SetLocalDerivative(3, fCosPhi * (fTrackPos[2] - fTrackPos0[2]));
765
766 // set global derivatives
3b2890be 767 SetGlobalDerivative(fDetElemNumber*fgNParCh+0, fSinPhi);
768 SetGlobalDerivative(fDetElemNumber*fgNParCh+1, -fCosPhi);
010eb601 769 if (fBFieldOn){
770 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
771 -fCosPhi*(fTrackPos[0]-fDetElemPos[0])
772 -fSinPhi*(fTrackPos[1]-fDetElemPos[1]));
773 }
774 else {
775 SetGlobalDerivative(fDetElemNumber*fgNParCh+2,
776 -fCosPhi*(fTrackPos0[0]+fTrackSlope0[0]*
777 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[0])
778 -fSinPhi*(fTrackPos0[1]+fTrackSlope0[1]*
779 (fTrackPos[2]-fTrackPos0[2])-fDetElemPos[1]));
780 }
781
782 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, fMeas[1], fSigma[1]);
783}
784
785void AliMUONAlignment::FillRecPointData() {
96ebe67e 786 /// Get information of current cluster
787 fClustPos[0] = fCluster->GetX();
788 fClustPos[1] = fCluster->GetY();
789 fClustPos[2] = fCluster->GetZ();
010eb601 790 fTransform->Global2Local(fDetElemId,fClustPos[0],fClustPos[1],fClustPos[2],
791 fClustPosLoc[0],fClustPosLoc[1],fClustPosLoc[2]);
792}
793
794void AliMUONAlignment::FillTrackParamData() {
96ebe67e 795 /// Get information of current track at current cluster
010eb601 796 fTrackPos[0] = fTrackParam->GetNonBendingCoor();
797 fTrackPos[1] = fTrackParam->GetBendingCoor();
798 fTrackPos[2] = fTrackParam->GetZ();
799 fTrackSlope[0] = fTrackParam->GetNonBendingSlope();
800 fTrackSlope[1] = fTrackParam->GetBendingSlope();
801 fTransform->Global2Local(fDetElemId,fTrackPos[0],fTrackPos[1],fTrackPos[2],
802 fTrackPosLoc[0],fTrackPosLoc[1],fTrackPosLoc[2]);
803}
804
805void AliMUONAlignment::FillDetElemData() {
806 /// Get information of current detection element
807 Double_t lDetElemLocX = 0.;
808 Double_t lDetElemLocY = 0.;
809 Double_t lDetElemLocZ = 0.;
96ebe67e 810 fDetElemId = fCluster->GetDetElemId();
010eb601 811 fDetElemNumber = fDetElemId%100;
812 for (int iCh=0; iCh<fDetElemId/100-1; iCh++){
813 fDetElemNumber += fgNDetElemCh[iCh];
814 }
815 fTransform->Local2Global(fDetElemId,lDetElemLocX,lDetElemLocY,lDetElemLocZ,
816 fDetElemPos[0],fDetElemPos[1],fDetElemPos[2]);
010eb601 817}
818
819void AliMUONAlignment::ProcessTrack(AliMUONTrack * track) {
96ebe67e 820 /// Process track; Loop over clusters and set local equations
010eb601 821 fTrack = track;
822 // get tclones arrays.
96ebe67e 823 fTrackParamAtCluster = fTrack->GetTrackParamAtCluster();
010eb601 824
825 // get size of arrays
96ebe67e 826 Int_t nTrackParam = fTrackParamAtCluster->GetEntries();
3b2890be 827 AliDebug(1,Form("Number of track param entries : %i ", nTrackParam));
010eb601 828
96ebe67e 829 for(Int_t iCluster=0; iCluster<nTrackParam; iCluster++) {
830 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtCluster()->At(iCluster);
831 fCluster = fTrackParam->GetClusterPtr();
832 if (!fCluster || !fTrackParam) continue;
010eb601 833 // fill local variables for this position --> one measurement
834 FillDetElemData();
835 FillRecPointData();
836 FillTrackParamData();
837// if (fDetElemId<500) continue;
838 fTrackPos0[0] = fTrackPos[0];
839 fTrackPos0[1] = fTrackPos[1];
840 fTrackPos0[2] = fTrackPos[2];
841 fTrackSlope0[0] = fTrackSlope[0];
842 fTrackSlope0[1] = fTrackSlope[1];
843 break;
844 }
845
96ebe67e 846 for(Int_t iCluster=0; iCluster<nTrackParam; iCluster++) {
010eb601 847 // and get new pointers
96ebe67e 848 fTrackParam = (AliMUONTrackParam *) fTrack->GetTrackParamAtCluster()->At(iCluster);
849 fCluster = fTrackParam->GetClusterPtr();
850 if (!fCluster || !fTrackParam) continue;
010eb601 851 // fill local variables for this position --> one measurement
852 FillDetElemData();
853 FillRecPointData();
854 FillTrackParamData();
855// if (fDetElemId<500) continue;
96ebe67e 856 AliDebug(1,Form("cluster: %i", iCluster));
010eb601 857 AliDebug(1,Form("x: %f\t y: %f\t z: %f\t DetElemID: %i\t ", fClustPos[0], fClustPos[1], fClustPos[2], fDetElemId));
858 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));
859
96ebe67e 860 AliDebug(1,Form("Track Parameter: %i", iCluster));
010eb601 861 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]));
862 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]));
863
864 fCosPhi = TMath::Cos(fPhi);
865 fSinPhi = TMath::Sin(fPhi);
866 if (fBFieldOn){
867 fMeas[0] = fTrackPos[0] - fClustPos[0];
868 fMeas[1] = fTrackPos[1] - fClustPos[1];
869 }
870 else {
871 fMeas[0] = - fClustPos[0];
872 fMeas[1] = - fClustPos[1];
873 }
874 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]));
875 // Set local equations
876 LocalEquationX();
877 LocalEquationY();
878 }
879}
880
881void AliMUONAlignment::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
882 /// Call local fit for this tracks
883 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
884 if (iRes && !lSingleFit) {
885 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
886 }
887}
888
889void AliMUONAlignment::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
890 /// Call global fit; Global parameters are stored in parameters
891 fMillepede->GlobalFit(parameters,errors,pulls);
3b2890be 892
010eb601 893 AliInfo("Done fitting global parameters!");
894 for (int iGlob=0; iGlob<fgNDetElem; iGlob++){
895 printf("%d\t %f\t %f\t %f \n",iGlob,parameters[iGlob*fgNParCh+0],parameters[iGlob*fgNParCh+1],parameters[iGlob*fgNParCh+2]);
896 }
897}
898
899Double_t AliMUONAlignment::GetParError(Int_t iPar) {
900 /// Get error of parameter iPar
901 Double_t lErr = fMillepede->GetParError(iPar);
902 return lErr;
903}
904
905void AliMUONAlignment::PrintGlobalParameters() {
906 /// Print global parameters
907 fMillepede->PrintGlobalParameters();
908}
909
910//_________________________________________________________________________
3b2890be 911TGeoCombiTrans AliMUONAlignment::ReAlign(const TGeoCombiTrans & transform, double *lMisAlignment) const
010eb601 912{
913 /// Realign given transformation by given misalignment and return the misaligned transformation
914
915 Double_t cartMisAlig[3] = {0,0,0};
916 Double_t angMisAlig[3] = {0,0,0};
3b2890be 917// const Double_t *trans = transform.GetTranslation();
918// TGeoRotation *rot;
919// // check if the rotation we obtain is not NULL
920// if (transform.GetRotation()) {
921// rot = transform.GetRotation();
922// }
923// else {
924// rot = new TGeoRotation("rot");
925// } // default constructor.
926
927 cartMisAlig[0] = -lMisAlignment[0];
928 cartMisAlig[1] = -lMisAlignment[1];
929 angMisAlig[2] = -lMisAlignment[2]*180./TMath::Pi();
930
931 TGeoTranslation deltaTrans(cartMisAlig[0], cartMisAlig[1], cartMisAlig[2]);
932 TGeoRotation deltaRot;
933 deltaRot.RotateX(angMisAlig[0]);
934 deltaRot.RotateY(angMisAlig[1]);
935 deltaRot.RotateZ(angMisAlig[2]);
936
937 TGeoCombiTrans deltaTransf(deltaTrans,deltaRot);
938 TGeoHMatrix newTransfMat = transform * deltaTransf;
939
940 return TGeoCombiTrans(newTransfMat);
010eb601 941}
942
943//______________________________________________________________________
944AliMUONGeometryTransformer *
945AliMUONAlignment::ReAlign(const AliMUONGeometryTransformer * transformer,
946 double *misAlignments, Bool_t verbose)
947
948{
949 /////////////////////////////////////////////////////////////////////
950 // Takes the internal geometry module transformers, copies them
951 // and gets the Detection Elements from them.
952 // Takes misalignment parameters and applies these
953 // to the local transform of the Detection Element
954 // Obtains the global transform by multiplying the module transformer
955 // transformation with the local transformation
956 // Applies the global transform to a new detection element
957 // Adds the new detection element to a new module transformer
958 // Adds the new module transformer to a new geometry transformer
959 // Returns the new geometry transformer
960
3b2890be 961 Double_t lModuleMisAlignment[3] = {0.,0.,0.};
010eb601 962 Double_t lDetElemMisAlignment[3] = {0.,0.,0.};
3b2890be 963 Int_t iDetElemId = 0;
964 Int_t iDetElemNumber = 0;
010eb601 965
966 AliMUONGeometryTransformer *newGeometryTransformer =
ef4cb4f1 967 new AliMUONGeometryTransformer();
3b2890be 968 for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
969 // module transformers
970 const AliMUONGeometryModuleTransformer *kModuleTransformer =
971 transformer->GetModuleTransformer(iMt, true);
010eb601 972
3b2890be 973 AliMUONGeometryModuleTransformer *newModuleTransformer =
974 new AliMUONGeometryModuleTransformer(iMt);
975 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
976
977 TGeoCombiTrans moduleTransform =
978 TGeoCombiTrans(*kModuleTransformer->GetTransformation());
979 // New module transformation
980 TGeoCombiTrans newModuleTransform = ReAlign(moduleTransform,lModuleMisAlignment);
981 newModuleTransformer->SetTransformation(newModuleTransform);
982
983 // Get delta transformation:
984 // Tdelta = Tnew * Told.inverse
985 TGeoHMatrix deltaModuleTransform =
986 AliMUONGeometryBuilder::Multiply(newModuleTransform,
987 kModuleTransformer->GetTransformation()->Inverse());
988 // Create module mis alignment matrix
989 newGeometryTransformer
990 ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
991
992 AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
993
994 if (verbose)
995 AliInfo(Form("%i DEs in old GeometryStore %i",detElements->GetSize(), iMt));
996
630711ed 997 TIter next(detElements->CreateIterator());
998 AliMUONGeometryDetElement* detElement;
999 Int_t iDe(-1);
1000 while ( ( detElement = static_cast<AliMUONGeometryDetElement*>(next()) ) )
1001 {
1002 ++iDe;
3b2890be 1003 /// make a new detection element
1004 AliMUONGeometryDetElement *newDetElement =
1005 new AliMUONGeometryDetElement(detElement->GetId(),
1006 detElement->GetVolumePath());
1007 TString lDetElemName(detElement->GetDEName());
1008 lDetElemName.ReplaceAll("DE","");
1009 iDetElemId = lDetElemName.Atoi();
1010 iDetElemNumber = iDetElemId%100;
1011 for (int iCh=0; iCh<iDetElemId/100-1; iCh++){
1012 iDetElemNumber += fgNDetElemCh[iCh];
1013 }
1014 for (int i=0; i<fgNParCh; i++) {
1015 lDetElemMisAlignment[i] = 0.0;
1016 if (iMt<fgNTrkMod) {
1017 AliInfo(Form("iMt %i, iCh %i, iDe %i, iDeId %i, iDeNb %i, iPar %i",iMt, iDetElemId/100, iDe, iDetElemId, iDetElemNumber, iDetElemNumber*fgNParCh+i));
1018 lDetElemMisAlignment[i] = misAlignments[iDetElemNumber*fgNParCh+i];
1019 }
1020 }
1021 // local transformation of this detection element.
1022 TGeoCombiTrans localTransform
1023 = TGeoCombiTrans(*detElement->GetLocalTransformation());
1024 TGeoCombiTrans newLocalTransform = ReAlign(localTransform,lDetElemMisAlignment);
1025 newDetElement->SetLocalTransformation(newLocalTransform);
1026
1027 // global transformation
1028 TGeoHMatrix newGlobalTransform =
1029 AliMUONGeometryBuilder::Multiply(newModuleTransform,
1030 newLocalTransform);
1031 newDetElement->SetGlobalTransformation(newGlobalTransform);
010eb601 1032
3b2890be 1033 // add this det element to module
1034 newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
1035 newDetElement);
1036
1037 // In the Alice Alignment Framework misalignment objects store
1038 // global delta transformation
1039 // Get detection "intermediate" global transformation
1040 TGeoHMatrix newOldGlobalTransform = newModuleTransform * localTransform;
1041 // Get detection element global delta transformation:
1042 // Tdelta = Tnew * Told.inverse
1043 TGeoHMatrix deltaGlobalTransform
1044 = AliMUONGeometryBuilder::Multiply(newGlobalTransform,
1045 newOldGlobalTransform.Inverse());
010eb601 1046
3b2890be 1047 // Create mis alignment matrix
1048 newGeometryTransformer
1049 ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
010eb601 1050 }
3b2890be 1051
1052 if (verbose)
1053 AliInfo(Form("Added module transformer %i to the transformer", iMt));
1054 newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
1055 }
010eb601 1056 return newGeometryTransformer;
1057}
1058
4818a9b7 1059//______________________________________________________________________
1060void AliMUONAlignment::SetAlignmentResolution(const TClonesArray* misAlignArray, Int_t rChId, Double_t rChResX, Double_t rChResY, Double_t rDeResX, Double_t rDeResY){
1061 //// Set alignment resolution to misalign objects to be stored in CDB
1062 Int_t chIdMin = (rChId<0)? 0 : rChId;
1063 Int_t chIdMax = (rChId<0)? 9 : rChId;
1064 Double_t chResX = rChResX;
1065 Double_t chResY = rChResY;
1066 Double_t deResX = rDeResX;
1067 Double_t deResY = rDeResY;
1068
1069 TMatrixDSym mChCorrMatrix(6);
1070 mChCorrMatrix[0][0]=chResX*chResX;
1071 mChCorrMatrix[1][1]=chResY*chResY;
1072 // mChCorrMatrix.Print();
1073
1074 TMatrixDSym mDECorrMatrix(6);
1075 mDECorrMatrix[0][0]=deResX*deResX;
1076 mDECorrMatrix[1][1]=deResY*deResY;
1077 // mDECorrMatrix.Print();
1078
1079 AliAlignObjMatrix *alignMat = 0x0;
1080
1081 for(Int_t chId=chIdMin; chId<=chIdMax; chId++) {
1082 TString chName1;
1083 TString chName2;
1084 if (chId<4){
1085 chName1 = Form("GM%d",chId);
1086 chName2 = Form("GM%d",chId);
1087 } else {
1088 chName1 = Form("GM%d",4+(chId-4)*2);
1089 chName2 = Form("GM%d",4+(chId-4)*2+1);
1090 }
1091
1092 for (int i=0; i<misAlignArray->GetEntries(); i++) {
1093 alignMat = (AliAlignObjMatrix*)misAlignArray->At(i);
1094 TString volName(alignMat->GetSymName());
1095 if((volName.Contains(chName1)&&
1096 ((volName.Last('/')==volName.Index(chName1)+chName1.Length())||
1097 (volName.Length()==volName.Index(chName1)+chName1.Length())))||
1098 (volName.Contains(chName2)&&
1099 ((volName.Last('/')==volName.Index(chName2)+chName2.Length())||
1100 (volName.Length()==volName.Index(chName2)+chName2.Length())))){
1101 volName.Remove(0,volName.Last('/')+1);
1102 if (volName.Contains("GM")) {
1103 // alignMat->Print("NULL");
1104 alignMat->SetCorrMatrix(mChCorrMatrix);
1105 } else if (volName.Contains("DE")) {
1106 // alignMat->Print("NULL");
1107 alignMat->SetCorrMatrix(mDECorrMatrix);
1108 }
1109 }
1110 }
1111 }
1112}