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