]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSAlignMille2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Implementing QA based on TPC online/offline tracks (Zhongbao)
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille2.cxx
CommitLineData
6be22b3f 1/**************************************************************************
2 * Copyright(c) 2007-2009, 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//-----------------------------------------------------------------------------
19//
20// Interface to AliMillePede2 alignment class for the ALICE ITS detector
21//
22// ITS specific alignment class which interface to AliMillepede.
23// For each track ProcessTrack calculates the local and global derivatives
24// at each hit and fill the corresponding local equations. Provide methods for
25// fixing or constraining detection elements for best results.
26//
27// author M. Lunardon (thanks to J. Castillo), ruben.shahoyan@cern.ch
28//-----------------------------------------------------------------------------
29
30#include <TFile.h>
31#include <TClonesArray.h>
32#include <TMath.h>
33#include <TVirtualFitter.h>
34#include <TGeoManager.h>
35#include <TSystem.h>
36#include <TRandom.h>
37#include <TCollection.h>
38#include <TGeoPhysicalNode.h>
39#include "AliITSAlignMille2.h"
40#include "AliITSgeomTGeo.h"
41#include "AliGeomManager.h"
42#include "AliMillePede2.h"
43#include "AliTrackPointArray.h"
44#include "AliAlignObjParams.h"
45#include "AliLog.h"
46#include "AliTrackFitterRieman.h"
47#include "AliITSAlignMille2Constraint.h"
48#include "AliITSAlignMille2ConstrArray.h"
49#include "AliITSresponseSDD.h"
50#include "AliITSTPArrayFit.h"
51#include "AliCDBManager.h"
52#include "AliCDBStorage.h"
53#include "AliCDBEntry.h"
54
55
56ClassImp(AliITSAlignMille2)
57
58const Char_t* AliITSAlignMille2::fgkRecKeys[] = {
59 "OCDB_PATH",
60 "GEOMETRY_FILE",
61 "SUPERMODULE_FILE",
62 "CONSTRAINTS_REFERENCE_FILE",
63 "PREALIGNMENT_FILE",
64 "PRECALIBSDD_FILE",
65 "INITCALBSDD_FILE",
66 "INITDELTA_FILE",
67 "SET_GLOBAL_DELTAS",
68 "CONSTRAINT_LOCAL",
69 "MODULE_VOLUID",
70 "MODULE_INDEX",
71 "SET_PSEUDO_PARENTS",
72 "SET_TRACK_FIT_METHOD",
73 "SET_MINPNT_TRA",
74 "SET_NSTDDEV",
75 "SET_RESCUT_INIT",
76 "SET_RESCUT_OTHER",
77 "SET_LOCALSIGMAFACTOR",
78 "SET_STARTFAC",
79 "SET_B_FIELD",
80 "SET_SPARSE_MATRIX",
81 "REQUIRE_POINT",
82 "CONSTRAINT_ORPHANS",
83 "CONSTRAINT_SUBUNITS",
84 "APPLY_CONSTRAINT",
85 "SET_EXTRA_CLUSTERS_MODE",
86 "SET_USE_TPAFITTER",
87 "SET_USE_LOCAL_YERROR",
88 "SET_MIN_POINTS_PER_MODULE"
89
90};
91
92const Char_t AliITSAlignMille2::fgkXYZ[] = "XYZ";
93
94//========================================================================================================
95
96AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;
97Int_t AliITSAlignMille2::fgInstanceID = 0;
98
99//________________________________________________________________________________________________________
100AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInfo )
101: TObject(),
102 fMillepede(0),
103 fStartFac(16.),
104 fResCutInitial(100.),
105 fResCut(100.),
106 fNGlobal(0),
107 fNLocal(4),
108 fNStdDev(3),
109 fIsMilleInit(kFALSE),
110 fAllowPseudoParents(kFALSE),
111 //
112 fTPAFitter(0),
113 fCurrentModule(0),
114 fTrack(0),
115 fTrackBuff(0),
116 fCluster(),
117 fCurrentSensID(-1),
118 fClusLoc(12*3),
119 fClusGlo(12*3),
120 fClusSigLoc(12*3),
121 fGlobalDerivatives(0),
122 fMeasLoc(0),
123 fMeasGlo(0),
124 fSigmaLoc(0),
125 fConstrPT(-1),
126 fConstrPTErr(-1),
127 fConstrCharge(0),
128 //
129 fMinNPtsPerTrack(3),
130 fInitTrackParamsMeth(1),
131 fTotBadLocEqPoints(0),
132 fRieman(0),
133 //
134 fConstraints(0),
135 fCacheMatrixOrig(kMaxITSSensID+1),
136 fCacheMatrixCurr(kMaxITSSensID+1),
137 //
138 fUseGlobalDelta(kFALSE),
139 fRequirePoints(kFALSE),
140 fTempExcludedModule(-1),
141 //
142 fDefCDBpath("local://$ALICE_ROOT/OCDB"),
143 fInitDeltaPath(""),
144 fInitSDDRespPath(""),
145 fPreCalSDDRespPath(""),
146 fGeometryPath("geometry.root"),
147 fPreDeltaPath(""),
148 fConstrRefPath(""),
149 fGeoManager(0),
150 fIsConfigured(kFALSE),
151 fPreAlignQF(0),
152//
153 fCorrectSDD(0),
154 fInitialRecSDD(0),
155 fPrealignment(0),
156 fConstrRef(0),
157 fMilleModule(2),
158 fSuperModule(2),
159 fNModules(0),
160 fNSuperModules(0),
161 fUsePreAlignment(kFALSE),
162 fUseLocalYErr(kFALSE),
163 fBOn(kFALSE),
164 fBField(0.0),
165 fMinPntPerSens(0),
166 fBug(0),
167 fMilleVersion(2),
168 fExtraClustersMode(0)
169{
170 /// main constructor that takes input from configuration file
171 for (int i=3;i--;) fSigmaFactor[i] = 1.0;
172 //
173 // new RS
174 for (Int_t i=0; i<6; i++) {
175 fNReqLayUp[i]=0;
176 fNReqLayDown[i]=0;
177 fNReqLay[i]=0;
178 }
179 for (Int_t i=0; i<3; i++) {
180 fNReqDetUp[i]=0;
181 fNReqDetDown[i]=0;
182 fNReqDet[i]=0;
183 }
184 //
185 if (ProcessUserInfo(userInfo)) exit(1);
186 //
187 Int_t lc=LoadConfig(configFilename);
188 if (lc) {
189 AliError(Form("Error %d loading configuration from %s",lc,configFilename));
190 exit(1);
191 }
192 //
193 fMillepede = new AliMillePede2();
194 fgInstance = this;
195 fgInstanceID++;
196 //
197}
198
199//________________________________________________________________________________________________________
200AliITSAlignMille2::~AliITSAlignMille2()
201{
202 /// Destructor
203 delete fMillepede;
204 delete[] fGlobalDerivatives;
205 delete fRieman;
206 delete fPrealignment;
207 delete fConstrRef;
208 delete fCorrectSDD;
209 delete fInitialRecSDD;
210 delete fTPAFitter;
211 fCacheMatrixOrig.Delete();
212 fCacheMatrixCurr.Delete();
213 fTrackBuff.Delete();
214 fConstraints.Delete();
215 fMilleModule.Delete();
216 fSuperModule.Delete();
217 if (--fgInstanceID==0) fgInstance = 0;
218}
219
220///////////////////////////////////////////////////////////////////////
221
222//________________________________________________________________________________________________________
223TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
224{
225 // read new record from config file
226 TString record;
227 static TObjArray* recElems = 0;
228 if (recElems) {delete recElems; recElems = 0;}
229 //
230 TString keyws = recTitle;
231 if (!keyws.IsNull()) {
232 keyws.ToUpper();
233 // keyws += " ";
234 }
235 while (record.Gets(stream)) {
236 int cmt=record.Index("#");
237 if (cmt>=0) record.Remove(cmt); // skip comment
238 record.ReplaceAll("\t"," ");
239 record.ReplaceAll("\r"," ");
240 record.Remove(TString::kBoth,' ');
241 if (record.IsNull()) continue; // nothing to decode
242 if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
243 //
244 recElems = record.Tokenize(" ");
245 recTitle = recElems->At(0)->GetName();
246 recTitle.ToUpper();
247 recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
248 break;
249 }
250 if (rew || !recElems) rewind(stream);
251 return recElems;
252}
253
254//________________________________________________________________________________________________________
255Int_t AliITSAlignMille2::CheckConfigRecords(FILE* stream)
256{
257 TString record,recTitle;
258 int lineCnt = 0;
259 rewind(stream);
260 while (record.Gets(stream)) {
261 int cmt=record.Index("#");
262 lineCnt++;
263 if (cmt>=0) record.Remove(cmt); // skip comment
264 record.ReplaceAll("\t"," ");
265 record.ReplaceAll("\r"," ");
266 record.Remove(TString::kBoth,' ');
267 if (record.IsNull()) continue; // nothing to decode
268 // extract keyword
269 int spc = record.Index(" ");
270 if (spc>0) recTitle = record(0,spc);
271 else recTitle = record;
272 recTitle.ToUpper();
273 Bool_t strOK = kFALSE;
274 for (int ik=kNKeyWords;ik--;) if (recTitle == fgkRecKeys[ik]) {strOK = kTRUE; break;}
275 if (strOK) continue;
276 //
277 AliError(Form("Unknown keyword %s at line %d",
278 recTitle.Data(),lineCnt));
279 return -1;
280 //
281 }
282 //
283 rewind(stream);
284 return 0;
285}
286
287
288//________________________________________________________________________________________________________
289Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
290{
291 // return 0 if success
292 // 1 if error in module index or voluid
293 //
294 FILE *pfc=fopen(cfile,"r");
295 if (!pfc) return -1;
296 //
297 TString record,recTitle,recOpt,recExt;
298 Int_t nrecElems,irec;
299 TObjArray *recArr=0;
300 //
301 fNModules = 0;
302 Bool_t stopped = kFALSE;
303 //
304 if (CheckConfigRecords(pfc)<0) return -1;
305 //
306 while(1) {
307 //
308 // ============= 1: we read some important records in predefined order ================
309 //
310 recTitle = fgkRecKeys[kOCDBPath];
311 if ( GetConfigRecord(pfc,recTitle,recOpt,1) && !(fDefCDBpath=recOpt).IsNull() ) {
312 AliInfo(Form("Configuration sets OCDB Def.Storage to %s",fDefCDBpath.Data()));
313 AliCDBManager::Instance()->SetDefaultStorage( fDefCDBpath.Data() );
314 }
315 //
316 //
317 recTitle = fgkRecKeys[kGeomFile];
318 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fGeometryPath = recOpt;
319 if ( InitGeometry() ) { AliError("Failed to find/load Geometry"); stopped = kTRUE; break;}
320 //
321 //
322 recTitle = fgkRecKeys[kSuperModileFile];
323 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
324 recOpt.IsNull() ||
325 gSystem->AccessPathName(recOpt.Data()) ||
326 LoadSuperModuleFile(recOpt.Data()))
327 { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
328 //
329 //
330 recTitle = fgkRecKeys[kConstrRefFile]; // LOCAL_CONSTRAINTS are defined wrt these deltas
331 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
332 if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
333 else {
334 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
335 if ( SetConstraintWrtRef(recOpt.Data()) )
336 { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
337 }
338 }
339 //
340 //
abd7ef79 341 recTitle = fgkRecKeys[kInitDeltaFile];
6be22b3f 342 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
343 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
abd7ef79 344 fInitDeltaPath = recOpt;
345 AliInfo(Form("Configuration sets Production Deltas to %s",fInitDeltaPath.Data()));
6be22b3f 346 }
6be22b3f 347 //
abd7ef79 348 // if initial deltas were provided, load them, apply to geometry and store are "original" matrices
349 if (CacheMatricesOrig()) {stopped = kTRUE; break;}
6be22b3f 350 //
abd7ef79 351 recTitle = fgkRecKeys[kPreDeltaFile];
6be22b3f 352 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
353 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
abd7ef79 354 fPreDeltaPath = recOpt;
355 AliInfo(Form("Configuration sets PreAlignment Deltas to %s",fPreDeltaPath.Data()));
6be22b3f 356 }
abd7ef79 357 if (LoadDeltas(fPreDeltaPath,fPrealignment)) {stopped = kTRUE; break;}
358 if (fPrealignment && ApplyToGeometry()) {stopped = kTRUE; break;}
6be22b3f 359 //
360 //
361 recTitle = fgkRecKeys[kPreCalSDDFile];
362 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull() ) {
363 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
364 fPreCalSDDRespPath = recOpt;
365 AliInfo(Form("Configuration sets PreCalibration SDD Response to %s",fPreCalSDDRespPath.Data()));
366 }
367 if (LoadSDDResponse(fPreCalSDDRespPath, fCorrectSDD) ) {stopped = kTRUE; break;}
368 //
369 recTitle = fgkRecKeys[ kInitCalSDDFile ];
370 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
371 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
372 fInitSDDRespPath = recOpt;
373 AliInfo(Form("Configuration sets Production SDD Response to %s",fInitSDDRespPath.Data()));
374 }
375 if (LoadSDDResponse(fInitSDDRespPath, fInitialRecSDD) ) {stopped = kTRUE; break;}
376 //
377 //
378 recTitle = fgkRecKeys[ kGlobalDeltas ];
379 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
380 //
381 // =========== 2: see if there are local gaussian constraints defined =====================
382 // Note that they should be loaded before the modules declaration
383 //
384 recTitle = fgkRecKeys[ kConstrLocal ];
385 while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
386 nrecElems = recArr->GetLast()+1;
387 if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
388 if (GetConstraint(recOpt.Data())) {
389 AliError(Form("Existing constraint %s repeated",recOpt.Data()));
390 stopped = kTRUE; break;
391 }
392 recExt = recArr->At(2)->GetName();
393 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
394 double val = recExt.Atof();
395 recExt = recArr->At(3)->GetName();
396 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
397 double err = recExt.Atof();
398 int nwgh = nrecElems - 4;
399 double *wgh = new double[nwgh];
400 for (nwgh=0,irec=4;irec<nrecElems;irec++) {
401 recExt = recArr->At(irec)->GetName();
402 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
403 wgh[nwgh++] = recExt.Atof();
404 }
405 if (stopped) {delete[] wgh; break;}
406 //
407 ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
408 delete[] wgh;
409 //
410 } // end while for loop over local constraints
411 if (stopped) break;
412 //
413 // =========== 3: now read modules to align ===================================
414 //
415 rewind(pfc);
416 while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
417 if (!(recTitle==fgkRecKeys[ kModVolID ] || recTitle==fgkRecKeys[ kModIndex ])) continue;
418 // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ] extra params]
419 // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
420 // sig* is the scaling parameters for the errors of the clusters of this module
421 // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
422 //
423 nrecElems = recArr->GetLast()+1;
424 if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
425 int idx = recOpt.Atoi();
426 UShort_t voluid = (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
427 AliITSAlignMille2Module* mod = 0;
428 //
429 if (voluid>=kMinITSSupeModuleID) { // custom supermodule
430 for (int j=0; j<fNSuperModules; j++) {
431 if (voluid==GetSuperModule(j)->GetVolumeID()) {
432 mod = new AliITSAlignMille2Module(*GetSuperModule(j));
433 // the matrix might be updated in case some prealignment was applied, check
434 TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
435 if (mup) *(mod->GetMatrix()) = *mup;
436 fMilleModule.AddAtAndExpand(mod,fNModules);
437 break;
438 }
439 }
440 }
441 else if (idx<=kMaxITSSensVID) {
442 mod = new AliITSAlignMille2Module(voluid);
443 fMilleModule.AddAtAndExpand(mod,fNModules);
444 }
445 if (!mod) {stopped = kTRUE; break;} // bad volid
446 //
447 // geometry variation settings
448 for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
449 irec = i+2;
450 if (irec >= nrecElems) break;
451 recExt = recArr->At(irec)->GetName();
452 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
453 mod->SetFreeDOF(i, recExt.Atof() );
454 }
455 if (stopped) break;
456 //
457 // scaling factors for cluster errors
458 // first set default ones
459 for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);
460 for (int i=0;i<3;i++) {
461 irec = i+8;
462 if (irec >= nrecElems) break;
463 recExt = recArr->At(irec)->GetName();
464 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
465 mod->SetSigmaFactor(i, recExt.Atof() );
466 }
467 if (stopped) break;
468 //
469 mod->SetGeomParamsGlobal(fUseGlobalDelta);
470 // now comes special detectors treatment
471 if (mod->IsSDD()) {
472 double vl = 0;
473 if (nrecElems>11) {
474 recExt = recArr->At(11)->GetName();
475 if (recExt.IsFloat()) vl = recExt.Atof();
476 else {stopped = kTRUE; break;}
477 irec = 11;
478 }
479 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
480 //
481 vl = 0;
482 if (nrecElems>12) {
483 recExt = recArr->At(12)->GetName();
484 if (recExt.IsFloat()) vl = recExt.Atof();
485 else {stopped = kTRUE; break;}
486 irec = 12;
487 }
488 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFDV,vl);
489 }
490 //
491 mod->SetUniqueID(fNModules);
492 mod->EvaluateDOF();
493 fNModules++;
494 //
495 // now check if there are local constraints on this module
496 for (++irec;irec<nrecElems;irec++) {
497 recExt = recArr->At(irec)->GetName();
498 if (recExt.IsFloat()) {stopped=kTRUE;break;}
499 AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
500 if (!cstr) {
501 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
502 stopped=kTRUE;
503 break;
504 }
505 cstr->AddModule(mod);
506 }
507 if (stopped) break;
508 } // end while for loop over modules
509 if (stopped) break;
510 //
511 if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}
512 BuildHierarchy(); // preprocess loaded modules
513 //
514 // =========== 4: the rest may come in arbitrary order =======================================
515 rewind(pfc);
516 while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
517 //
518 nrecElems = recArr->GetLast()+1;
519 //
520 // some simple flags -----------------------------------------------------------------------
521 //
522 if (recTitle == fgkRecKeys[ kPseudoParents ]) SetAllowPseudoParents(kTRUE);
523 //
524 // some optional parameters ----------------------------------------------------------------
525 else if (recTitle == fgkRecKeys[ kTrackFitMethod ]) {
526 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
527 SetInitTrackParamsMeth(recOpt.Atoi());
528 }
529 //
530 else if (recTitle == fgkRecKeys[ kMinPntTrack ]) {
531 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
532 fMinNPtsPerTrack = recOpt.Atoi();
533 }
534 //
535 else if (recTitle == fgkRecKeys[ kNStDev ]) {
536 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
537 fNStdDev = (Int_t)recOpt.Atof();
538 }
539 //
540 else if (recTitle == fgkRecKeys[ kResCutInit ]) {
541 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
542 fResCutInitial = recOpt.Atof();
543 }
544 //
545 else if (recTitle == fgkRecKeys[ kResCutOther ]) {
546 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
547 fResCut = recOpt.Atof();
548 }
549 //
550 else if (recTitle == fgkRecKeys[ kLocalSigFactor ]) { //-------------------------
551 for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
552 fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
553 if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
554 }
555 if (stopped) break;
556 }
557 //
558 else if (recTitle == fgkRecKeys[ kStartFactor ]) { //-------------------------
559 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
560 fStartFac = recOpt.Atof();
561 }
562 //
563 // pepo2708909
564 else if (recTitle == fgkRecKeys[ kExtraClustersMode ]) { //-------------------------
565 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
566 fExtraClustersMode = recOpt.Atoi();
567 }
568 // endpepo270809
569 //
570 else if (recTitle == fgkRecKeys[ kBField ]) { //-------------------------
571 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
572 SetBField( recOpt.Atof() );
573 }
574 //
575 else if (recTitle == fgkRecKeys[ kSparseMatrix ]) { // matrix solver type
576 //
577 AliMillePede2::SetGlobalMatSparse(kTRUE);
578 if (recOpt.IsNull()) continue;
579 // solver type and settings
580 if (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
581 else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
582 else {stopped = kTRUE; break;}
583 //
584 if (nrecElems>=3) { // preconditioner type
585 recExt = recArr->At(2)->GetName();
586 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
587 AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
588 }
589 //
590 if (nrecElems>=4) { // tolerance
591 recExt = recArr->At(3)->GetName();
592 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
593 AliMillePede2::SetMinResTol( recExt.Atof() );
594 }
595 //
596 if (nrecElems>=5) { // maxIter
597 recExt = recArr->At(4)->GetName();
598 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
599 AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
600 }
601 }
602 //
603 else if (recTitle == fgkRecKeys[ kRequirePoint ]) { //-------------------------
604 // syntax: REQUIRE_POINT where ndet updw nreqpts
605 // where = LAYER or DETECTOR
606 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
607 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
608 // nreqpts = minimum number of points of that type
609 if (nrecElems>=5) {
610 recOpt.ToUpper();
611 int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
612 int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
613 int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
614 fRequirePoints = kTRUE;
615 if (recOpt == "LAYER") {
616 if (lr<0 || lr>5) {stopped = kTRUE; break;}
617 if (hb>0) fNReqLayUp[lr] = np;
618 else if (hb<0) fNReqLayDown[lr] = np;
619 else fNReqLay[lr] = np;
620 }
621 else if (recOpt == "DETECTOR") {
622 if (lr<0 || lr>2) {stopped = kTRUE; break;}
623 if (hb>0) fNReqDetUp[lr] = np;
624 else if (hb<0) fNReqDetDown[lr] = np;
625 else fNReqDet[lr] = np;
626 }
627 else {stopped = kTRUE; break;}
628 }
629 else {stopped = kTRUE; break;}
630 }
631 //
632 // global constraints on the subunits/orphans
633 else if (recTitle == fgkRecKeys[ kConstrOrphans ]) { //------------------------
634 // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
635 if (nrecElems<4) {stopped = kTRUE; break;}
636 recExt = recArr->At(2)->GetName();
637 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
638 double val = recExt.Atof();
639 UInt_t pattern = 0;
640 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
641 recExt = recArr->At(irec)->GetName();
642 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
643 pattern |= 0x1 << recExt.Atoi();
644 }
645 if (stopped) break;
646 if (recOpt == "MEAN") ConstrainOrphansMean(val,pattern);
647 else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
648 else {stopped = kTRUE; break;}
649 }
650 //
651 else if (recTitle == fgkRecKeys[ kConstrSubunits ]) { //------------------------
652 // expect ONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
653 if (nrecElems<5) {stopped = kTRUE; break;}
654 recExt = recArr->At(2)->GetName();
655 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
656 double val = recExt.Atof();
657 UInt_t pattern = 0;
658 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
659 recExt = recArr->At(irec)->GetName();
660 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
661 int parid = recExt.Atoi();
662 if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
663 else break; // list of params is over
664 }
665 if (stopped) break;
666 //
667 Bool_t meanC;
668 if (recOpt == "MEAN") meanC = kTRUE;
669 else if (recOpt == "MEDIAN") meanC = kFALSE;
670 else {stopped = kTRUE; break;}
671 //
672 int curID = -1;
673 int rangeStart = -1;
674 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
675 recExt = recArr->At(irec)->GetName();
676 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
677 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
678 else curID = recExt.Atoi();
679 //
680 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
681 // this was a range start or single
682 int start;
683 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
684 else start = curID; // create constraint either for single module (or 1st in the range)
685 for (int id=start;id<=curID;id++) {
686 int id0 = IsVIDDefined(id);
687 if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
688 if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
689 else ConstrainModuleSubUnitsMedian(id0,val,pattern);
690 }
691 }
692 if (rangeStart>=0) stopped = kTRUE; // unfinished range
693 if (stopped) break;
694 }
695 //
696 // association of modules with local constraints
697 else if (recTitle == fgkRecKeys[ kApplyConstr ]) { //------------------------
698 // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
699 if (nrecElems<3) {stopped = kTRUE; break;}
700 int nmID0=-1,nmID1=-1;
701 for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
702 recExt = recArr->At(irec)->GetName();
703 if (recExt.IsFloat()) break;
704 // check if such a constraint was declared
705 if (!GetConstraint(recExt.Data())) {
706 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
707 stopped=kTRUE;
708 break;
709 }
710 if (nmID0<0) nmID0 = irec;
711 nmID1 = irec;
712 }
713 if (stopped) break;
714 //
715 if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
716 //
717 // now read the list of modules to constrain
718 int curID = -1;
719 int rangeStart = -1;
720 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
721 recExt = recArr->At(irec)->GetName();
722 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
723 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
724 else curID = recExt.Atoi();
725 //
726 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
727 //
728 // this was a range start or single
729 int start;
730 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
731 else start = curID; // create constraint either for single module (or 1st in the range)
732 for (int id=start;id<=curID;id++) {
733 AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
734 if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
735 for (int nmid=nmID0;nmid<=nmID1;nmid++)
736 ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
737 }
738 }
739 if (rangeStart>=0) stopped = kTRUE; // unfinished range
740 if (stopped) break;
741 }
742 // Do we use new TrackPointArray fitter ?
743 else if (recTitle == fgkRecKeys[ kTPAFitter ]) {
744 // expect SET_TPAFITTER
745 fTPAFitter = new AliITSTPArrayFit(kNLocal);
746 }
747 // Do we use new local Y errors?
748 else if (recTitle == fgkRecKeys[ kUseLocalYErr ]) {
749 // expect SET_TPAFITTER
750 fUseLocalYErr = kTRUE;
751 }
752 //
753 else if (recTitle == fgkRecKeys[ kMinPointsSens ]) { //-------------------------
754 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
755 SetMinPointsPerSensor( recOpt.Atoi() );
756 }
757 //
758 else continue; // already processed record
759 //
760 } // end of while loop 4 over the various params
761 //
762 break;
763 } // end of while(1) loop
764 //
765 fclose(pfc);
766 if (stopped) {
767 AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
768 return -1;
769 }
770 //
abd7ef79 771 if (CacheMatricesCurr()) return -1;
6be22b3f 772 SetUseLocalYErrors(fUseLocalYErr); // YErr used only with TPAFitter
773 fIsConfigured = kTRUE;
774 return 0;
775}
776
777//________________________________________________________________________________________________________
778void AliITSAlignMille2::BuildHierarchy()
779{
780 // build the hieararhy of the modules to align
781 //
782 if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
783 AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
784 "Since Deltas are local, switching to NoPseudoParents");
785 SetAllowPseudoParents(kFALSE);
786 }
787 // set parent/child relationship for modules to align
788 AliInfo("Setting parent/child relationships\n");
789 //
790 // 1) child -> parent reference
791 for (int ipar=0;ipar<fNModules;ipar++) {
792 AliITSAlignMille2Module* parent = GetMilleModule(ipar);
793 if (parent->IsSensor()) continue; // sensor cannot be a parent
794 //
795 for (int icld=0;icld<fNModules;icld++) {
796 if (icld==ipar) continue;
797 AliITSAlignMille2Module* child = GetMilleModule(icld);
798 if (!child->BelongsTo(parent)) continue;
799 // child cannot have more sensors than the parent
800 if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
801 //
802 AliITSAlignMille2Module* parOld = child->GetParent();
803 // is this parent candidate closer than the old parent ?
804 if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
805 child->SetParent(parent);
806 }
807 //
808 }
809 //
810 // add parent -> children reference
811 for (int icld=0;icld<fNModules;icld++) {
812 AliITSAlignMille2Module* child = GetMilleModule(icld);
813 AliITSAlignMille2Module* parent = child->GetParent();
814 if (parent) parent->AddChild(child);
815 }
816 //
817 // reorder the modules in such a way that parents come first
818 for (int icld=0;icld<fNModules;icld++) {
819 AliITSAlignMille2Module* child = GetMilleModule(icld);
820 AliITSAlignMille2Module* parent;
821 while ( (parent=child->GetParent()) && (parent->GetUniqueID()>child->GetUniqueID()) ) {
822 // swap
823 fMilleModule[icld] = parent;
824 fMilleModule[parent->GetUniqueID()] = child;
825 child->SetUniqueID(parent->GetUniqueID());
826 parent->SetUniqueID(icld);
827 child = parent;
828 }
829 //
830 }
831 //
832 // Go over the child->parent chain and mark modules with explicitly provided sensors.
833 // If the sensors of the unit are explicitly declared, all undeclared sensors are
834 // suppresed in this unit.
835 for (int icld=fNModules;icld--;) {
836 AliITSAlignMille2Module* child = GetMilleModule(icld);
837 AliITSAlignMille2Module* parent = child->GetParent();
838 if (!parent) continue;
839 //
840 // check if this parent was already processed
841 if (!parent->AreSensorsProvided()) {
842 parent->DelSensitiveVolumes();
843 parent->SetSensorsProvided(kTRUE);
844 }
845 // reattach sensors to parent
846 for (int isc=child->GetNSensitiveVolumes();isc--;) {
847 UShort_t senVID = child->GetSensVolVolumeID(isc);
848 if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
849 }
850 }
851 //
852}
853
854// pepo
855//________________________________________________________________________________________________________
856void AliITSAlignMille2::SetCurrentModule(Int_t id)
857{
858 // set the current supermodule
859 // new meaning
860 if (fMilleVersion>=2) {
861 fCurrentModule = GetMilleModule(id);
862 return;
863 }
864 // old meaning
865 if (fMilleVersion<=1) {
866 Int_t index=id;
867 /// set as current the SuperModule that contains the 'index' sens.vol.
868 if (index<0 || index>2197) {
869 AliInfo("index does not correspond to a sensitive volume!");
870 return;
871 }
872 UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
873 Int_t k=IsContained(voluid);
874 if (k>=0){
875 fCurrentSensID = index;
876 fCluster.SetVolumeID(voluid);
877 fCluster.SetXYZ(0,0,0);
878 InitModuleParams();
879 }
880 else
881 AliInfo(Form("module %d not defined\n",index));
882 }
883}
884// endpepo
885//________________________________________________________________________________________________________
886void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts)
887{
888 // set minimum number of points in specific detector or layer
889 // where = LAYER or DETECTOR
890 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
891 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
892 // nreqpts = minimum number of points of that type
893 ndet--;
894 if (strstr(where,"LAYER")) {
895 if (ndet<0 || ndet>5) return;
896 if (updw>0) fNReqLayUp[ndet]=nreqpts;
897 else if (updw<0) fNReqLayDown[ndet]=nreqpts;
898 else fNReqLay[ndet]=nreqpts;
899 fRequirePoints=kTRUE;
900 }
901 else if (strstr(where,"DETECTOR")) {
902 if (ndet<0 || ndet>2) return;
903 if (updw>0) fNReqDetUp[ndet]=nreqpts;
904 else if (updw<0) fNReqDetDown[ndet]=nreqpts;
905 else fNReqDet[ndet]=nreqpts;
906 fRequirePoints=kTRUE;
907 }
908}
909
910//________________________________________________________________________________________________________
911Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname)
912{
913 /// index from symname
914 if (!symname) return -1;
915 for (Int_t i=0;i<=kMaxITSSensID; i++) {
916 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
917 }
918 return -1;
919}
920
921//________________________________________________________________________________________________________
922Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid)
923{
924 /// index from volume ID
925 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
926 if (lay<1|| lay>6) return -1;
927 Int_t idx=Int_t(voluid)-2048*lay;
928 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
929 for (Int_t ilay=1; ilay<lay; ilay++)
930 idx += AliGeomManager::LayerSize(ilay);
931 return idx;
932}
933
934//________________________________________________________________________________________________________
935UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname)
936{
937 /// volume ID from symname
938 /// works for sensitive volumes only
939 if (!symname) return 0;
940
941 for (UShort_t voluid=2000; voluid<13300; voluid++) {
942 Int_t modId;
943 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
944 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
945 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
946 }
947 }
948
949 return 0;
950}
951
952//________________________________________________________________________________________________________
953UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index)
954{
955 /// volume ID from index
956 if (index<0) return 0;
957 if (index<2198)
958 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
959 else {
960 for (int i=0; i<fNSuperModules; i++) {
961 if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
962 }
963 }
964 return 0;
965}
966
967//________________________________________________________________________________________________________
968Int_t AliITSAlignMille2::InitGeometry()
969{
970 /// initialize geometry
971 AliInfo("Loading initial geometry");
972 if (!fGeometryPath.IsNull() && gSystem->AccessPathName(fGeometryPath.Data()) ) {
973 AliError(Form("Explicitly provided geometry file %s is not accessible",fGeometryPath.Data()));
974 return -1;
975 }
976 //
977 AliGeomManager::LoadGeometry(fGeometryPath.Data());
978 fGeoManager = AliGeomManager::GetGeometry();
979 if (!fGeoManager) {
980 AliInfo("Couldn't initialize geometry");
981 return -1;
982 }
983 return 0;
984}
985
986//________________________________________________________________________________________________________
987Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname)
988{
989 // Load the global deltas from this file. The local gaussian constraints on some modules
990 // will be defined with respect to the deltas from this reference file, converted to local
991 // delta format. Note: conversion to local format requires reloading the geometry!
992 //
993 AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
994 if (!fGeoManager) return -1;
995 fConstrRefPath = reffname;
996 if (fConstrRefPath == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
997 fConstrRef = new TClonesArray("AliAlignObjParams",1);
998 return 0;
999 }
1000 if (LoadDeltas(fConstrRefPath,fConstrRef)) return -1;
1001 //
1002 // we need ideal geometry to convert global deltas to local ones
1003 if (fUsePreAlignment) {
1004 AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
1005 return -1;
1006 }
1007 //
1008 AliInfo("Converting global reference deltas to local ones");
1009 Int_t nprea = fConstrRef->GetEntriesFast();
1010 for (int ix=0; ix<nprea; ix++) {
1011 AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
1012 if (!preo->ApplyToGeometry()) return -1;
1013 }
1014 //
1015 // now convert the global reference deltas to local ones
1016 for (int i=fConstrRef->GetEntriesFast();i--;) {
1017 AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
1018 TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
1019 if (!mupd) { // this is not alignable entry, need to look in the supermodules
1020 for (int im=fNSuperModules;im--;) {
1021 AliITSAlignMille2Module* mod = GetSuperModule(im);
1022 if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
1023 mupd = mod->GetMatrix();
1024 break;
1025 }
1026 if (!mupd) {
1027 AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
1028 return -1;
1029 }
1030 }
1031 TGeoHMatrix preMat;
1032 preo->GetMatrix(preMat); // Delta_Glob
1033 TGeoHMatrix tmpMat = *mupd; // Delta_Glob * Delta_Glob_Par * M
1034 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
1035 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
1036 preo->SetMatrix(tmpMat); // local corrections
1037 }
1038 //
1039 // we need to reload the geometry spoiled by this reference deltas...
1040 delete fGeoManager;
1041 AliInfo("Reloading initial geometry");
1042 return InitGeometry();
1043 //
1044}
1045
1046//________________________________________________________________________________________________________
1047void AliITSAlignMille2::Init()
1048{
1049 // perform global initialization
1050 //
1051 if (fIsMilleInit) {
1052 AliInfo("Millepede has been already initialized!");
1053 return;
1054 }
1055 // range constraints in such a way that the childs are constrained before their parents
1056 // orphan constraints come last
1057 for (int ic=0;ic<GetNConstraints();ic++) {
1058 for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
1059 AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
1060 AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
1061 if (cst0->GetModuleID()<cst1->GetModuleID()) {
1062 // swap
1063 fConstraints[ic] = cst1;
1064 fConstraints[ic1] = cst0;
1065 }
1066 }
1067 }
1068 //
1069 if (!GetUseGlobalDelta()) {
1070 AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
1071 for (int imd=fNModules;imd--;) {
1072 AliITSAlignMille2Module* mod = GetMilleModule(imd);
1073 int npar = mod->GetNParTot();
1074 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1075 for (int ipar=0;ipar<npar;ipar++) {
1076 if (!mod->IsFreeDOF(ipar)) continue;
1077 mod->SetParOffset(ipar,fNGlobal++);
1078 }
1079 }
1080 }
1081 else {
1082 // init millepede, decide which parameters are to be fitted explicitly
1083 for (int imd=fNModules;imd--;) {
1084 AliITSAlignMille2Module* mod = GetMilleModule(imd);
1085 int npar = mod->GetNParTot();
1086 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1087 for (int ipar=0;ipar<npar;ipar++) {
1088 if (!mod->IsFreeDOF(ipar)) continue; // fixed
1089 //
1090 int nFreeInstances = 0;
1091 //
1092 AliITSAlignMille2Module* parent = mod;
1093 Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
1094 //
1095 Bool_t addToFit = kFALSE;
1096 // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
1097 // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
1098 // 2) the same applies to all of its parents
1099 // 3) it has at least 1 unconstrained direct child
1100 while(parent) {
1101 if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
1102 nFreeInstances++;
1103 if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
1104 if (cstGauss) addToFit = kTRUE;
1105 parent = parent->GetParent();
1106 }
1107 if (nFreeInstances>1) {
1108 AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
1109 "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
1110 exit(1);
1111 }
1112 //
1113 // i) Are PseudoParents allowed?
1114 if (!PseudoParentsAllowed()) addToFit = kTRUE;
1115 // ii) check if this module has no child with such a free parameter. Since the order of this check
1116 // goes from child to parent, by this moment such a parameter must have been already added
1117 else if (!IsParModFamilyVaried(mod,ipar)) addToFit = kTRUE; // no varied children at all
1118 else if (!IsParFamilyFree(mod,ipar,1)) addToFit = kTRUE; // no unconstrained direct children
1119 // otherwise the value of this parameter can be extracted from simple contraint and the values of
1120 // the relevant parameters of its children the fit is done. Hence it is not included
1121 if (!addToFit) continue;
1122 //
1123 // shall add this parameter to explicit fit
1124 // printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
1125 mod->SetParOffset(ipar,fNGlobal++);
1126 }
1127 }
1128 }
1129 //
1130 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, kNLocal, fNStdDev));
1131 fGlobalDerivatives = new Double_t[fNGlobal];
1132 memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1133 //
1134 fMillepede->InitMille(fNGlobal,kNLocal,fNStdDev,fResCut,fResCutInitial);
1135 fMillepede->SetMinPntValid(fMinPntPerSens);
1136 fIsMilleInit = kTRUE;
1137 //
1138 ResetLocalEquation();
1139 AliInfo("Parameters initialized to zero");
1140 //
1141 /// Fix non free parameters
1142 for (Int_t i=0; i<fNModules; i++) {
1143 AliITSAlignMille2Module* mod = GetMilleModule(i);
1144 for (Int_t j=0; j<mod->GetNParTot(); j++) {
1145 if (mod->GetParOffset(j)<0) continue; // not varied
1146 FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1147 fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1148 }
1149 }
1150 //
1151 // Set iterations
1152 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
1153 //
1154}
1155
1156//________________________________________________________________________________________________________
1157void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma)
1158{
1159 /// Constrain equation defined by par to value
1160 if (!fIsMilleInit) Init();
1161 fMillepede->SetGlobalConstraint(par, value, sigma);
1162 AliInfo("Adding constraint");
1163}
1164
1165//________________________________________________________________________________________________________
1166void AliITSAlignMille2::InitGlobalParameters(Double_t *par)
1167{
1168 /// Initialize global parameters with par array
1169 if (!fIsMilleInit) Init();
1170 fMillepede->SetGlobalParameters(par);
1171 AliInfo("Init Global Parameters");
1172}
1173
1174//________________________________________________________________________________________________________
1175void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value)
1176{
1177 /// Parameter iPar is encourage to vary in [-value;value].
1178 /// If value == 0, parameter is fixed
1179 if (!fIsMilleInit) {
1180 AliInfo("Millepede has not been initialized!");
1181 return;
1182 }
1183 fMillepede->SetParSigma(iPar, value);
1184 if (IsZero(value)) AliInfo(Form("Parameter %i Fixed", iPar));
1185}
1186
1187//________________________________________________________________________________________________________
1188void AliITSAlignMille2::ResetLocalEquation()
1189{
1190 /// Reset the derivative vectors
1191 for(int i=kNLocal;i--;) fLocalDerivatives[i] = 0.0;
1192 memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
1193}
1194
1195//________________________________________________________________________________________________________
1196Int_t AliITSAlignMille2::ApplyToGeometry()
1197{
1198 // apply prealignment to ideal geometry
1199 Int_t nprea = fPrealignment->GetEntriesFast();
1200 AliInfo(Form("Array of prealignment deltas: %d entries",nprea));
1201 //
1202 for (int ix=0; ix<nprea; ix++) {
1203 AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1204 Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1205 if (index>=0) {
1206 if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1207 fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1208 }
1209 if (!preo->ApplyToGeometry()) {
1210 AliError(Form("Failed on ApplyToGeometry at %s",preo->GetSymName()));
1211 return -1;
1212 }
1213 }
1214 //
1215 fUsePreAlignment = kTRUE;
1216 return 0;
1217}
1218
1219//________________________________________________________________________________________________________
1220Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1221{
1222 // quality factors from prealignment
1223 if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1224 return fPreAlignQF[index]-1;
1225}
1226
1227//________________________________________________________________________________________________________
1228AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp)
1229{
1230 /// create a new AliTrackPointArray keeping only defined modules
1231 /// move points according to a given prealignment, if any
1232 /// sort alitrackpoints w.r.t. global Y direction, if selected
1233 const Double_t kRad2L[6] = {5*5,10*10,18*18,30*30,40*40,60*60};
1234 const Float_t kSensSigY2[6] = {200e-4*200e-4/12, 200e-4*200e-4/12,
1235 300e-4*300e-4/12, 300e-4*300e-4/12,
1236 300e-4*300e-4/12, 300e-4*300e-4/12}; // thickness^2/12
1237 //
1238 fTrack = NULL;
1239 Int_t idx[20];
1240 Short_t lrID[20];
1241 Int_t npts=atp->GetNPoints();
1242 TGeoHMatrix hcov;
1243 //
1244 /// checks if AliTrackPoints belong to defined modules
1245 Int_t ngoodpts=0;
1246 Int_t intidx[20];
1247 for (int j=0; j<npts; j++) {
1248 intidx[j] = IsVIDContained(atp->GetVolumeID()[j]);
1249 if (intidx[j]<0) continue;
1250 ngoodpts++;
1251 Float_t xx=atp->GetX()[j];
1252 Float_t yy=atp->GetY()[j];
1253 Float_t r=xx*xx + yy*yy;
1254 int lay;
1255 for (lay=0;lay<6;lay++) if (r<kRad2L[lay]) break;
1256 if (lay>5) continue;
1257 lrID[j] = lay;
1258 }
1259 //
1260 AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
1261
1262 // pepo270809
1263 Int_t nextra=0;
1264 // extra clusters selection mode
1265 if (fExtraClustersMode) {
1266 // 1 = keep one cluster, remove randomly the extra
1267 // 2 = keep one cluster, remove the internal one
1268 // 10 = keep tracks only if at least one extra is present
1269
1270 int iextra1[20],iextra2[20],layovl[20];
1271 // extra clusters mapping
1272 for (Int_t ipt=0; ipt<npts; ipt++) {
1273 if (intidx[ipt]<0) continue; // looks only defined modules...
1274 float p1x=atp->GetX()[ipt];
1275 float p1y=atp->GetY()[ipt];
1276 float p1z=atp->GetZ()[ipt];
1277 int lay1=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ipt]));
1278 float r1 = p1x*p1x + p1y*p1y;
1279 UShort_t volid1=atp->GetVolumeID()[ipt];
1280
1281 for (int ik=ipt+1; ik<npts; ik++) {
1282 if (intidx[ik]<0) continue;
1283 // compare point ipt with next ones
1284 int lay2=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ik]));
1285 // check if same layer
1286 if (lay2 != lay1) continue;
1287 UShort_t volid2=atp->GetVolumeID()[ik];
1288 // check if different module
1289 if (volid1 == volid2) continue;
1290
1291 float p2x=atp->GetX()[ik];
1292 float p2y=atp->GetY()[ik];
1293 float p2z=atp->GetZ()[ik];
1294 float r2 = p2x*p2x + p2y*p2y;
1295 float dr= (p1x-p2x)*(p1x-p2x) + (p1y-p2y)*(p1y-p2y) + (p1z-p2z)*(p1z-p2z);
1296
1297 // looks for pairs with dr<1 cm, same layer but different module
1298 if (dr<1.0) {
1299 // extra1 is the one with smaller radius in rphi plane
1300 if (r1<r2) {
1301 iextra1[nextra]=ipt;
1302 iextra2[nextra]=ik;
1303 }
1304 else {
1305 iextra1[nextra]=ik;
1306 iextra2[nextra]=ipt;
1307 }
1308 layovl[nextra]=lay1;
1309 nextra++;
1310 }
1311 }
1312 } // end overlaps mapping
1313
1314 // mode=1: keep only one clusters and remove the other randomly
1315 if (fExtraClustersMode==1 && nextra) {
1316 for (int ie=0; ie<nextra; ie++) {
1317 if (gRandom->Rndm()<0.5)
1318 intidx[iextra1[ie]]=-1;
1319 else
1320 intidx[iextra2[ie]]=-1;
1321 }
1322 }
1323
1324 // mode=2: keep only one clusters and remove the other...
1325 if (fExtraClustersMode==2 && nextra) {
1326 for (int ie=0; ie<nextra; ie++) {
1327 if (layovl[ie]==1) intidx[iextra2[ie]]=-1;
1328 else if (layovl[ie]==2) intidx[iextra1[ie]]=-1;
1329 else intidx[iextra1[ie]]=-1;
1330 }
1331 }
1332
1333 // mode=10: reject track if no overlaps are present
1334 if (fExtraClustersMode==10 && nextra==0) {
1335 AliInfo("Track with no extra clusters: rejected!");
1336 return NULL;
1337 }
1338
1339 // recalculate ngoodpts
1340 ngoodpts=0;
1341 for (int i=0; i<npts; i++) {
1342 if (intidx[i]>=0) ngoodpts++;
1343 }
1344 }
1345 // endpepo270809
1346
1347 // reject track if not enough points are left
1348 if (ngoodpts<fMinNPtsPerTrack) {
1349 AliInfo("Track with not enough points!");
1350 return NULL;
1351 }
1352 // >> RS
1353 AliTrackPoint p;
1354 // check points in specific places
1355 if (fRequirePoints) {
1356 Int_t nlayup[6],nlaydown[6],nlay[6];
1357 Int_t ndetup[3],ndetdown[3],ndet[3];
1358 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1359 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1360
1361 for (int i=0; i<npts; i++) {
1362 // skip not defined points
1363 if (intidx[i]<0) continue;
1364 //
1365 Float_t yy=atp->GetY()[i];
1366 int lay = lrID[i];
1367 int det=lay/2;
1368 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
1369
1370 if (yy>=0.0) { // UP point
1371 nlayup[lay]++;
1372 nlay[lay]++;
1373 ndetup[det]++;
1374 ndet[det]++;
1375 }
1376 else {
1377 nlaydown[lay]++;
1378 nlay[lay]++;
1379 ndetdown[det]++;
1380 ndet[det]++;
1381 }
1382 }
1383 //
1384 // checks minimum values
1385 Bool_t isok=kTRUE;
1386 for (Int_t j=0; j<6; j++) {
1387 if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE;
1388 if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE;
1389 if (nlay[j]<fNReqLay[j]) isok=kFALSE;
1390 }
1391 for (Int_t j=0; j<3; j++) {
1392 if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE;
1393 if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE;
1394 if (ndet[j]<fNReqDet[j]) isok=kFALSE;
1395 }
1396 if (!isok) {
1397 AliDebug(2,Form("Track does not meet all location point requirements!"));
1398 return NULL;
1399 }
1400 }
1401 // build a new track with (sorted) (prealigned) good points
1402 // pepo200709
1403 //fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
1404 fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts];
1405 if (!fTrack) {
1406 fTrack = new AliTrackPointArray(ngoodpts);
1407 // fTrackBuff.AddAtAndExpand(fTrack,ngoodpts-fMinNPtsPerTrack);
1408 fTrackBuff.AddAtAndExpand(fTrack,ngoodpts);
1409 }
1410 // fTrack = new AliTrackPointArray(ngoodpts);
1411 // endpepo200709
1412 //
1413 //
1414 for (int i=0; i<npts; i++) idx[i]=i;
1415 // sort track if required
1416 TMath::Sort(npts,atp->GetY(),idx); // sort descending...
1417 //
1418 Int_t npto=0;
1419 if (fClusLoc.GetSize()<3*npts) fClusLoc.Set(3*npts);
1420 if (fClusGlo.GetSize()<3*npts) fClusGlo.Set(3*npts);
1421 if (fClusSigLoc.GetSize()<3*npts) fClusSigLoc.Set(3*npts);
1422 //
1423 for (int i=0; i<npts; i++) {
1424 // skip not defined points
1425 if (intidx[idx[i]]<0) continue;
1426 atp->GetPoint(p,idx[i]);
1427 int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1428 //
1429 // prealign point if required
1430 // get matrix used to produce the digits
1431 AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1432 TGeoHMatrix *svOrigMatrix = GetSensorOrigMatrixSID(sid); //mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1433 // get back real local coordinate
1434 Double_t *pl = fClusLoc.GetArray() + npto*3;
1435 Double_t *pg = fClusGlo.GetArray() + npto*3;
1436 Double_t *sgl = fClusSigLoc.GetArray() + npto*3;
1437 pg[0]=p.GetX();
1438 pg[1]=p.GetY();
1439 pg[2]=p.GetZ();
1440 AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1441 svOrigMatrix->MasterToLocal(pg,pl);
abd7ef79 1442 AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]));
6be22b3f 1443 //
1444 // this is a temporary code to extract the drift speed used for given point
1445 if (p.GetDriftTime()>0) { // RRR
1446 // calculate the drift speed
1447 fDriftTime0[npto] = fInitialRecSDD ? fInitialRecSDD->GetTimeZero(sid) : 0.;
1448 double tdif = p.GetDriftTime() - fDriftTime0[npto];
1449 if (tdif<=0) tdif = 1;
1450 double vdrift = (3.5085-TMath::Abs(pl[0]))/tdif;
1451 if (vdrift<0) vdrift = 0;
1452 //
1453 // TEMPORARY CORRECTION (if provided) -------------->>>
1454 if (fCorrectSDD) {
1455 float t0Upd = fCorrectSDD->GetTimeZero(sid);
1456 vdrift += fCorrectSDD->GetDeltaVDrift(sid);
1457 tdif = p.GetDriftTime() - t0Upd;
1458 // correct Xlocal
1459 pl[0] = TMath::Sign(3.5085 - vdrift*tdif,pl[0]);
1460 fDriftTime0[npto] = t0Upd;
1461 }
1462 // TEMPORARY CORRECTION (if provided) --------------<<<
1463 fDriftSpeed[npto] = TMath::Sign(vdrift,pl[0]);
1464 //
1465 }
1466 // update covariance matrix
1467 Double_t hcovel[9];
1468 hcovel[0]=double(p.GetCov()[0]);
1469 hcovel[1]=double(p.GetCov()[1]);
1470 hcovel[2]=double(p.GetCov()[2]);
1471 hcovel[3]=double(p.GetCov()[1]);
1472 hcovel[4]=double(p.GetCov()[3]);
1473 hcovel[5]=double(p.GetCov()[4]);
1474 hcovel[6]=double(p.GetCov()[2]);
1475 hcovel[7]=double(p.GetCov()[4]);
1476 hcovel[8]=double(p.GetCov()[5]);
1477 hcov.SetRotation(hcovel);
abd7ef79 1478 //
1479 if (AliLog::GetGlobalDebugLevel()>=2) {
1480 AliInfo("Original Global Cov Matrix");
1481 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovel[0],hcovel[1],hcovel[2],hcovel[4],hcovel[5],hcovel[8]);
1482 }
1483 //
6be22b3f 1484 // now rotate in local system
1485 hcov.Multiply(svOrigMatrix);
1486 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1487 // now hcov is LOCAL COVARIANCE MATRIX
1488 // apply sigma scaling
abd7ef79 1489 Double_t *hcovscl = hcov.GetRotationMatrix();
1490 if (AliLog::GetGlobalDebugLevel()>=2) {
1491 AliInfo("Original Local Cov Matrix");
1492 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1493 }
1494 hcovscl[4] = fUseLocalYErr ? kSensSigY2[lrID[idx[i]]] : 1E-8; // error due to the sensor thickness
6be22b3f 1495 //
1496 for (int ir=3;ir--;) for (int ic=3;ic--;) {
1497 if (ir==ic) {
abd7ef79 1498 if ( IsZero(hcovscl[ir*3+ic],1e-8) ) hcovscl[ir*3+ic] = 1E-8;
6be22b3f 1499 else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1500 sgl[ir] = TMath::Sqrt(hcovscl[ir*3+ic]);
1501 }
1502 else hcovscl[ir*3+ic] = 0;
1503 }
1504 //
abd7ef79 1505 if (AliLog::GetGlobalDebugLevel()>=2) {
1506 AliInfo("Modified Local Cov Matrix");
1507 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1508 }
6be22b3f 1509 //
1510 if (fBug==1) {
1511 // correzione bug LAYER 5 SSD temporanea..
1512 int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1513 if (ssdidx>=500 && ssdidx<1248) {
1514 int ladder=(ssdidx-500)%22;
1515 if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1516 if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1517 }
1518 }
1519 /// get (evenctually prealigned) matrix of sens. vol.
1520 TGeoHMatrix *svMatrix = GetSensorCurrMatrixSID(sid); //mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
1521 // modify global coordinates according with pre-aligment
1522 svMatrix->LocalToMaster(pl,pg);
1523 // now rotate in local system
1524 hcov.Multiply(&svMatrix->Inverse());
1525 hcov.MultiplyLeft(svMatrix); // hcov is back in GLOBAL RF
1526 // cure once more
1527 for (int ir=3;ir--;) for (int ic=3;ic--;) if (IsZero(hcovscl[ir*3+ic])) hcovscl[ir*3+ic] = 0.;
1528 // printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1529 //
abd7ef79 1530 if (AliLog::GetGlobalDebugLevel()>=2) {
1531 AliInfo("Modified Global Cov Matrix");
1532 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1533 }
1534 //
6be22b3f 1535 Float_t pcov[6];
1536 pcov[0]=hcovscl[0];
1537 pcov[1]=hcovscl[1];
1538 pcov[2]=hcovscl[2];
1539 pcov[3]=hcovscl[4];
1540 pcov[4]=hcovscl[5];
1541 pcov[5]=hcovscl[8];
1542 //
1543 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1544 // printf("New Gl coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1545 AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1546 fTrack->AddPoint(npto,&p);
1547 AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,fTrack->GetX()[npto],
1548 fTrack->GetY()[npto],fTrack->GetZ()[npto],fTrack->GetVolumeID()[npto] ));
1549 // printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY());
1550 npto++;
1551 }
1552 //
1553 return fTrack;
1554}
1555
1556//________________________________________________________________________________________________________
1557AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp)
1558{
1559 /// sort alitrackpoints w.r.t. global Y direction
1560 AliTrackPointArray *atps=NULL;
1561 Int_t idx[20];
1562 Int_t npts=atp->GetNPoints();
1563 AliTrackPoint p;
1564 atps=new AliTrackPointArray(npts);
1565
1566 TMath::Sort(npts,atp->GetY(),idx);
1567
1568 for (int i=0; i<npts; i++) {
1569 atp->GetPoint(p,idx[i]);
1570 atps->AddPoint(i,&p);
1571 AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
1572 }
1573 return atps;
1574}
1575
1576//________________________________________________________________________________________________________
1577Int_t AliITSAlignMille2::GetCurrentLayer() const
1578{
1579 // get current layer id
1580 if (!fGeoManager) {
1581 AliInfo("ITS geometry not initialized!");
1582 return -1;
1583 }
1584 return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1585}
1586
1587//________________________________________________________________________________________________________
1588Int_t AliITSAlignMille2::InitModuleParams()
1589{
1590 /// initialize geometry parameters for a given detector
1591 /// for current cluster (fCluster)
1592 /// fGlobalInitParam[] is set as:
1593 /// [tx,ty,tz,psi,theta,phi]
1594 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1595 /// *** At the moment: using Raffalele's angles definition ***
1596 ///
1597 /// return 0 if success
1598 /// If module is found but has no parameters to vary, return 1
1599
1600 if (!fGeoManager) {
1601 AliInfo("ITS geometry not initialized!");
1602 return -1;
1603 }
1604
1605 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1606
1607 // set the internal index (index in module list)
1608 UShort_t voluid=fCluster.GetVolumeID();
1609 fCurrentSensID = AliITSAlignMille2Module::GetIndexFromVolumeID(voluid);
1610 //
1611 // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1612 Int_t k=fNModules-1;
1613 fCurrentModule = 0;
1614 // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules
1615 while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
1616 if (k<0) return -3;
1617 //
1618 for (int i=AliITSAlignMille2Module::kMaxParTot;i--;) fModuleInitParam[i] = 0.0;
1619 //
1620 int clID = fCluster.GetUniqueID()-1;
1621 if (clID<0) { // external cluster
1622 fMeasGlo = &fExtClusterPar[0];
1623 fMeasLoc = &fExtClusterPar[3];
1624 fSigmaLoc = &fExtClusterPar[6];
1625 fExtClusterPar[0] = fCluster.GetX();
1626 fExtClusterPar[1] = fCluster.GetY();
1627 fExtClusterPar[2] = fCluster.GetZ();
1628 //
1629 TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1630 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1631 TGeoHMatrix hcov;
1632 Double_t hcovel[9];
1633 hcovel[0]=double(fCluster.GetCov()[0]);
1634 hcovel[1]=double(fCluster.GetCov()[1]);
1635 hcovel[2]=double(fCluster.GetCov()[2]);
1636 hcovel[3]=double(fCluster.GetCov()[1]);
1637 hcovel[4]=double(fCluster.GetCov()[3]);
1638 hcovel[5]=double(fCluster.GetCov()[4]);
1639 hcovel[6]=double(fCluster.GetCov()[2]);
1640 hcovel[7]=double(fCluster.GetCov()[4]);
1641 hcovel[8]=double(fCluster.GetCov()[5]);
1642 hcov.SetRotation(hcovel);
1643 // now rotate in local system
1644 hcov.Multiply(svMatrix);
1645 hcov.MultiplyLeft(&svMatrix->Inverse());
1646 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
1647 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
1648 //
1649 }
1650 else {
1651 int offs = 3*clID;
1652 fMeasGlo = fClusGlo.GetArray() + offs;
1653 fMeasLoc = fClusLoc.GetArray() + offs;
1654 fSigmaLoc = fClusSigLoc.GetArray() + offs;
1655 }
1656 //
1657 // set minimum value for SigmaLoc to 10 micron
1658 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
1659 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
1660 //
1661 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
1662 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
1663 //
1664 return 0;
1665}
1666
1667//________________________________________________________________________________________________________
1668void AliITSAlignMille2::Print(Option_t*) const
1669{
1670 // print current status
1671 printf("*** AliMillepede for ITS ***\n");
1672 printf(" Number of defined super modules: %d\n",fNModules);
1673 printf(" Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
1674 //
1675 if (fGeoManager)
1676 printf(" geometry loaded from %s\n",fGeometryPath.Data());
1677 else
1678 printf(" geometry not loaded\n");
1679 //
1680 if (fUsePreAlignment)
1681 printf(" using prealignment from %s \n",fPreDeltaPath.Data());
1682 else
1683 printf(" prealignment not used\n");
1684 //
1685 //
1686 if (fBOn)
1687 printf(" B Field set to %f T - using helices\n",fBField);
1688 else
1689 printf(" B Field OFF - using straight lines \n");
1690 //
1691 if (fTPAFitter)
1692 printf(" Using AliITSTPArrayFit class for track fitting\n");
1693 else
1694 printf(" Using StraightLine/Riemann fitter for track fitting\n");
1695 //
1696 printf("Using local Y error due to the sensor thickness: %s\n",(fUseLocalYErr && fTPAFitter) ? "ON":"OFF");
1697 //
1698 if (fRequirePoints) printf(" Required points in tracks:\n");
1699 for (Int_t i=0; i<6; i++) {
1700 if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
1701 if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
1702 if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]);
1703 }
1704 for (Int_t i=0; i<3; i++) {
1705 if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
1706 if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
1707 if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]);
1708 }
1709 //
1710 printf("\n Millepede configuration parameters:\n");
1711 printf(" init value for chi2 cut : %.4f\n",fStartFac);
1712 printf(" first iteration cut value : %.4f\n",fResCutInitial);
1713 printf(" other iterations cut value : %.4f\n",fResCut);
1714 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
1715 printf(" def.scaling for local sigmas : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
1716 printf(" min.tracks per module : %d\n",fMinPntPerSens);
1717 //
1718 printf("List of defined modules:\n");
1719 printf(" intidx\tindex\tvoluid\tname\n");
1720 for (int i=0; i<fNModules; i++) {
1721 AliITSAlignMille2Module* md = GetMilleModule(i);
1722 printf(" %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
1723 }
1724}
1725
1726//________________________________________________________________________________________________________
1727AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
1728{
1729 // return pointer to a defined supermodule
1730 // return NULL if error
1731 Int_t i=IsVIDDefined(voluid);
1732 if (i<0) return NULL;
1733 return GetMilleModule(i);
1734}
1735
1736//________________________________________________________________________________________________________
1737AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
1738{
1739 // return pointer to a defined supermodule
1740 // return NULL if error
1741 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1742 if (vid>0) return GetMilleModuleByVID(vid);
1743 else { // this is not alignable module, need to look within defined supermodules
1744 int i = IsSymDefined(symname);
1745 if (i>=0) return GetMilleModule(i);
1746 }
1747 return 0;
1748}
1749
1750//________________________________________________________________________________________________________
1751AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
1752{
1753 // return pointer to a defined/contained supermodule
1754 // return NULL otherwise
1755 int i = IsSymContained(symname);
1756 return i<0 ? 0 : GetMilleModule(i);
1757}
1758
1759//________________________________________________________________________________________________________
1760AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
1761{
1762 // get delta from prealignment for given volume
1763 if (!fPrealignment) return 0;
1764 for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1765 AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
1766 if (!strcmp(preob->GetSymName(),symname)) return preob;
1767 }
1768 return 0;
1769}
1770
1771//________________________________________________________________________________________________________
1772AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
1773{
1774 // get delta with respect to which the constraint is declared
1775 if (!fConstrRef) return 0;
1776 for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1777 AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
1778 if (!strcmp(preob->GetSymName(),symname)) return preob;
1779 }
1780 return 0;
1781}
1782
1783//________________________________________________________________________________________________________
1784Bool_t AliITSAlignMille2::InitRiemanFit()
1785{
1786 // Initialize Riemann Fitter for current track
1787 // return kFALSE if error
1788
1789 if (!fBOn) return kFALSE;
1790
1791 Int_t npts=0;
1792 AliTrackPoint ap;
1793 npts = fTrack->GetNPoints();
1794 AliDebug(3,Form("Fitting track with %d points",npts));
1795 if (!fRieman) fRieman = new AliTrackFitterRieman();
1796 fRieman->Reset();
1797 fRieman->SetTrackPointArray(fTrack);
1798
1799 TArrayI ai(npts);
1800 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
1801
1802 // fit track with 5 params in his own tracking-rotated reference system
1803 // xc = -p[1]/p[0];
1804 // yc = 1/p[0];
1805 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
1806 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
1807 return kFALSE;
1808 }
1809
1810 for (int i=0; i<5; i++)
1811 fLocalInitParam[i] = fRieman->GetParam()[i];
1812
1813 return kTRUE;
1814}
1815
1816//________________________________________________________________________________________________________
1817void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int flag)
1818{
1819 // local function for minuit
1820 const double kTiny = 1.e-14;
1821 chi2 = 0;
1822 static AliTrackPoint pnt;
1823 static Bool_t fullErr2D;
1824 //
1825 if (flag==1) fullErr2D = kFALSE;//kTRUE;
1826 fullErr2D = kTRUE;
1827 enum {kAX,kAZ,kBX,kBZ};
1828 enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1829 //
1830 AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
1831 AliTrackPointArray* track = alig->GetCurrentTrack();
1832 //
1833 int npts = track->GetNPoints();
1834 for (int ip=0;ip<npts;ip++) {
1835 track->GetPoint(pnt,ip);
1836 const float *cov = pnt.GetCov();
1837 double y = pnt.GetY();
1838 double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
1839 double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
1840 double xxe = cov[kXX];
1841 double zze = cov[kZZ];
1842 double xze = cov[kXZ];
1843 //
1844 if (fullErr2D) {
1845 xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
1846 zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
1847 xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
1848 }
1849 //
1850 double det = xxe*zze - xze*xze;
1851 if (det<kTiny) {
1852 printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
1853 "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
1854 xxe = cov[kXX];
1855 zze = cov[kZZ];
1856 xze = cov[kXZ];
1857 fullErr2D = kFALSE;
1858 }
1859 double xxeI = zze/det;
1860 double zzeI = xxe/det;
1861 double xzeI =-xze/det;
1862 //
1863 chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
1864 //
1865 // printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI, chi2);
1866 }
1867 //
1868}
1869
1870//________________________________________________________________________________________________________
1871void AliITSAlignMille2::InitTrackParams(int meth)
1872{
1873 /// initialize local parameters with different methods
1874 /// for current track (fTrack)
1875 Int_t npts=0;
1876 AliTrackPoint ap;
1877 double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
1878 // simple linear interpolation
1879 // get local starting parameters (to be substituted by ESD track parms)
1880 // local parms (fLocalInitParam[]) are:
1881 // [0] = global x coord. of straight line intersection at y=0 plane
1882 // [1] = global z coord. of straight line intersection at y=0 plane
1883 // [2] = px/py
1884 // [3] = pz/py
1885 // test #1: linear fit in x(y) and z(y)
1886 npts = fTrack->GetNPoints();
1887 AliDebug(3,Form("*** initializing track with %d points ***",npts));
1888 for (int i=npts;i--;) {
1889 sY += fTrack->GetY()[i];
1890 sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
1891 sX += fTrack->GetX()[i];
1892 sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
1893 sZ += fTrack->GetZ()[i];
1894 sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
1895 }
1896 det = sYY*npts-sY*sY;
1897 if (IsZero(det)) det = 1E-16;
1898 fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
1899 fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
1900 //
1901 fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
1902 fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
1903 // pepo200709
1904 fLocalInitParam[4] = 0.0;
1905 // endpepo200709
1906
1907 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f ugx = %f\n",fLocalInitParam[0],fLocalInitParam[2]));
1908 //
1909 if (meth==1) return;
1910 //
1911 // perform full fit accounting for cov.matrix
1912 static TVirtualFitter *minuit = 0;
1913 static Double_t step[5] = {1E-3,1E-3,1E-4,1E-4,1E-5};
1914 static Double_t arglist[10];
1915 //
1916 if (!minuit) {
1917 minuit = TVirtualFitter::Fitter(0,4);
1918 minuit->SetFCN(trackFit2D);
1919 arglist[0] = 1;
1920 minuit->ExecuteCommand("SET ERR",arglist, 1);
1921 //
1922 arglist[0] = -1;
1923 minuit->ExecuteCommand("SET PRINT",arglist,1);
1924 //
1925 }
1926 //
1927 minuit->SetParameter(0, "ax", fLocalInitParam[0], step[0], 0,0);
1928 minuit->SetParameter(1, "az", fLocalInitParam[1], step[1], 0,0);
1929 minuit->SetParameter(2, "bx", fLocalInitParam[2], step[2], 0,0);
1930 minuit->SetParameter(3, "bz", fLocalInitParam[3], step[3], 0,0);
1931 //
1932 arglist[0] = 1000; // number of function calls
1933 arglist[1] = 0.001; // tolerance
1934 minuit->ExecuteCommand("MIGRAD",arglist,2);
1935 //
1936 for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
1937 for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
1938 /*
1939 double amin,edm,errdef;
1940 int nvpar,nparx;
1941 minuit->GetStats(amin,edm,errdef,nvpar,nparx);
1942 amin /= (2*npts - 4);
1943 printf("Mchi2: %+e\n",amin);
1944 */
1945 //
1946}
1947
1948//________________________________________________________________________________________________________
1949Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
1950{
1951 // checks if supermodule with this symname is defined and return the internal index
1952 // return -1 if not.
1953 for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
1954 return -1;
1955}
1956
1957//________________________________________________________________________________________________________
1958Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
1959{
1960 // checks if module with this symname is defined and return the internal index
1961 // return -1 if not.
1962 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1963 if (vid>0) return IsVIDContained(vid);
1964 // only sensors have real vid, but maybe we have a supermodule with fake vid?
1965 // IMPORTANT: always start from the end to start from the sensors
1966 return IsSymDefined(symname);
1967}
1968
1969//________________________________________________________________________________________________________
1970Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
1971{
1972 // checks if supermodule 'voluid' is defined and return the internal index
1973 // return -1 if not.
1974 for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
1975 return -1;
1976}
1977
1978//________________________________________________________________________________________________________
1979Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
1980{
1981 // checks if the sensitive module 'voluid' is contained inside a supermodule
1982 // and return the internal index of the last identified supermodule
1983 // return -1 if error
1984 // IMPORTANT: always start from the end to start from the sensors
1985 if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
1986 for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
1987 return -1;
1988}
1989
1990//________________________________________________________________________________________________________
1991Int_t AliITSAlignMille2::CheckCurrentTrack()
1992{
1993 /// checks if AliTrackPoints belongs to defined modules
1994 /// return number of good poins
1995 /// return 0 if not enough points
1996
1997 Int_t npts = fTrack->GetNPoints();
1998 Int_t ngoodpts=0;
1999 // debug points
2000 for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
2001 //
2002 if (ngoodpts<fMinNPtsPerTrack) return 0;
2003
2004 return ngoodpts;
2005}
2006
2007//________________________________________________________________________________________________________
2008Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track)
2009{
2010 /// Process track; Loop over hits and set local equations
2011 /// here 'track' is a AliTrackPointArray
2012 /// return 0 if success;
2013
2014 if (!fIsMilleInit) Init();
2015 //
2016 Int_t npts = track->GetNPoints();
2017 AliDebug(2,Form("*** Input track with %d points ***",npts));
2018
2019 // preprocessing of the input track: keep only points in defined volumes,
2020 // move points if prealignment is set, sort by Yglo if required
2021
2022 fTrack=PrepareTrack(track);
2023 if (!fTrack) return -1;
2024
2025 npts = fTrack->GetNPoints();
2026 if (npts>kMaxPoints) {
2027 AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
2028 }
2029 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
2030 //
2031 if (fTPAFitter) { // use dediacted fitter
2032 //
2033 fTPAFitter->AttachPoints(fTrack);
2034 if (fBOn) fTPAFitter->SetBz(fBField);
2035 if (fInitTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
2036 double chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2037 //
2038 // suppress eventual constraints to not affect fit of the next track
2039 fConstrCharge = 0;
2040 fConstrPT = fConstrPTErr = -1;
2041 //
2042 if ( chi2<0 || (fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) {
2043 AliInfo("Track fit failed! skipping this track...");
2044 fTPAFitter->Reset();
2045 fTrack = NULL;
2046 return -5;
2047 }
2048 /*
2049 double *pr = fTPAFitter->GetParams();
2050 printf("FtPar: %+.5e %+.5e %+.5e %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
2051 */
2052 }
2053 else {
2054 //
2055 if (!fBOn) { // straight lines
2056 // set local starting parameters (to be substituted by ESD track parms)
2057 // local parms (fLocalInitParam[]) are:
2058 // [0] = global x coord. of straight line intersection at y=0 plane
2059 // [1] = global z coord. of straight line intersection at y=0 plane
2060 // [2] = px/py
2061 // [3] = pz/py
2062 InitTrackParams(fInitTrackParamsMeth);
2063 /*
2064 double *pr = fLocalInitParam;
2065 printf("FtPar: %+.5e %+.5e %+.5e %+.5e |\n",pr[0],pr[1],pr[2],pr[3]); // RRR
2066 */
2067 }
2068 else {
2069 // local parms (fLocalInitParam[]) are the Riemann Fitter params
2070 if (!InitRiemanFit()) {
2071 AliInfo("Riemann fit failed! skipping this track...");
2072 fTrack=NULL;
2073 return -5;
2074 }
2075 }
2076 }
2077 //
2078 // printf("Params: "); for (int i=0;i<fNLocal;i++) printf("%+.2e ",fLocalInitParam[i]); printf("\n");//RRR
2079 Int_t nloceq=0;
2080 Int_t ngloeq=0;
2081 static Mille2Data md[kMaxPoints];
2082 //
2083 for (Int_t ipt=0; ipt<npts; ipt++) {
2084 fTrack->GetPoint(fCluster,ipt);
2085 fCluster.SetUniqueID(ipt+1);
2086 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
2087
2088 // set geometry parameters for the the current module
2089 if (InitModuleParams()) continue;
2090 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n",
2091 track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
2092 fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
2093 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
2094 int res = fTPAFitter ? AddLocalEquationTPA(md[nloceq]) : AddLocalEquation(md[nloceq]);
2095 if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
2096 else if (res>0) {nloceq++; ngloeq++;}
2097 } // end loop over points
2098 //
2099 fTrack=NULL;
2100 // not enough good points?
2101 if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
2102 //
2103 // finally send local equations to millepede
2104 SetLocalEquations(md,nloceq);
2105 fMillepede->SaveRecordData(); // RRR
2106 //
2107 return 0;
2108}
2109
2110//________________________________________________________________________________________________________
2111Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar)
2112{
2113 /// calculate track intersection point in local coordinates
2114 /// according with a given set of parameters (local(4) and global(6))
2115 /// and fill fPintLoc/Glo
2116 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
2117 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
2118 /// return 0 if success
2119
2120 AliDebug(3,Form("lpar = %g %g %g %g %g\ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],lpar[4],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
2121 AliDebug(3,Form("deltalpar = %g %g %g %g %g\n",lpar[0]-fLocalInitParam[0],lpar[1]-fLocalInitParam[1],lpar[2]-fLocalInitParam[2],lpar[3]-fLocalInitParam[3],lpar[4]-fLocalInitParam[4]));
2122
2123
2124 // prepare the TGeoHMatrix
2125 TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
2126 !fUseGlobalDelta);
2127 if (!tempHMat) return -1;
2128
2129 Double_t v0g[3]; // vector with straight line direction in global coord.
2130 Double_t p0g[3]; // point of the straight line (glo)
2131
2132 if (fBOn) { // B FIELD!
2133 AliTrackPoint prf;
2134 for (int ip=0; ip<5; ip++)
2135 fRieman->SetParam(ip,lpar[ip]);
2136
2137 if (!fRieman->GetPCA(fCluster,prf)) {
2138 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
2139 return -3;
2140 }
2141 // now determine straight line passing tangent to fit curve at prf
2142 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
2143 // mo' P1=(X,Y,Z)_glo_prf
2144 // => (x,y,Z)_trk_prf ruotando di alpha...
2145 Double_t alpha=fRieman->GetAlpha();
2146 Double_t x1g = prf.GetX();
2147 Double_t y1g = prf.GetY();
2148 Double_t z1g = prf.GetZ();
2149 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
2150 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
2151 Double_t z1t = z1g;
2152
2153 Double_t x2t = x1t+1.0;
2154 Double_t y2t = y1t+fRieman->GetDYat(x1t);
2155 Double_t z2t = z1t+fRieman->GetDZat(x1t);
2156 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
2157 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
2158 Double_t z2g = z2t;
2159
2160 AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi()));
2161 AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
2162 AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
2163
2164 if (TMath::Abs(y2g-y1g)<1e-15) {
2165 AliInfo("DeltaY=0! Cannot proceed...");
2166 return -1;
2167 }
2168 // ugx,1,ugz
2169 v0g[0] = (x2g-x1g)/(y2g-y1g);
2170 v0g[1]=1.0;
2171 v0g[2] = (z2g-z1g)/(y2g-y1g);
2172
2173 // point: just keep prf
2174 p0g[0]=x1g;
2175 p0g[1]=y1g;
2176 p0g[2]=z1g;
2177 }
2178 else { // staight line
2179 // vector of initial straight line direction in glob. coord
2180 v0g[0]=lpar[2];
2181 v0g[1]=1.0;
2182 v0g[2]=lpar[3];
2183
2184 // intercept in yg=0 plane in glob coord
2185 p0g[0]=lpar[0];
2186 p0g[1]=0.0;
2187 p0g[2]=lpar[1];
2188 }
2189 AliDebug(3,Form("Line vector: ( %f , %f , %f ) point:( %f , %f , %f )\n",v0g[0],v0g[1],v0g[2],p0g[0],p0g[1],p0g[2]));
2190
2191 // same in local coord.
2192 Double_t p0l[3],v0l[3];
2193 tempHMat->MasterToLocalVect(v0g,v0l);
2194 tempHMat->MasterToLocal(p0g,p0l);
2195
2196 if (TMath::Abs(v0l[1])<1e-15) {
2197 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2198 return -1;
2199 }
2200
2201 // local intersection point
2202 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2203 fPintLoc[1] = 0;
2204 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2205
2206 // global intersection point
2207 tempHMat->LocalToMaster(fPintLoc,fPintGlo);
2208 AliDebug(3,Form("Intesect. point: L( %f , %f , %f ) G( %f , %f , %f )\n",fPintLoc[0],fPintLoc[1],fPintLoc[2],fPintGlo[0],fPintGlo[1],fPintGlo[2]));
2209
2210 return 0;
2211}
2212
2213//________________________________________________________________________________________________________
2214Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
2215{
2216 /// calculate numerically (ROOT's style) the derivatives for
2217 /// local X intersection and local Z intersection
2218 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
2219 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2220 /// return 0 if success
2221
2222 // copy initial parameters
2223 Double_t lpar[kNLocal];
2224 Double_t gpar[kNParCh];
2225 Double_t *derivative;
2226 for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2227 for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
2228
2229 // trial with fixed dpar...
2230 Double_t dpar = 0.;
2231
2232 if (islpar) { // track parameters
2233 //dpar=fLocalInitParam[paridx]*0.001;
2234 // set minimum dpar
2235 derivative = fDerivativeLoc[paridx];
2236 if (!fBOn) {
2237 if (paridx<3) dpar=1.0e-4; // translations
2238 else dpar=1.0e-6; // direction
2239 }
2240 else { // B Field
2241 // pepo: proviamo con 1/1000, poi evenctually 1/100...
2242 Double_t dfrac=0.01;
2243 switch(paridx) {
2244 case 0:
2245 // RMS cosmics: 1e-4
2246 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2247 break;
2248 case 1:
2249 // RMS cosmics: 0.2
2250 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2251 break;
2252 case 2:
2253 // RMS cosmics: 9
2254 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2255 break;
2256 case 3:
2257 // RMS cosmics: 7
2258 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2259 break;
2260 case 4:
2261 // RMS cosmics: 0.3
2262 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2263 break;
2264 }
2265 }
2266 }
2267 else { // alignment global parameters
2268 derivative = fDerivativeGlo[paridx];
2269 //dpar=fModuleInitParam[paridx]*0.001;
2270 if (paridx<3) dpar=1.0e-4; // translations
2271 else dpar=1.0e-2; // angles
2272 }
2273
2274 AliDebug(3,Form("+++ using dpar=%g",dpar));
2275
2276 // calculate derivative ROOT's like:
2277 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2278 Double_t pintl1[3]; // f(x-h)
2279 Double_t pintl2[3]; // f(x-h/2)
2280 Double_t pintl3[3]; // f(x+h/2)
2281 Double_t pintl4[3]; // f(x+h)
2282
2283 // first values
2284 if (islpar) lpar[paridx] -= dpar;
2285 else gpar[paridx] -= dpar;
2286 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2287 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2288
2289 // second values
2290 if (islpar) lpar[paridx] += dpar/2;
2291 else gpar[paridx] += dpar/2;
2292 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2293 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2294
2295 // third values
2296 if (islpar) lpar[paridx] += dpar;
2297 else gpar[paridx] += dpar;
2298 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2299 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2300
2301 // fourth values
2302 if (islpar) lpar[paridx] += dpar/2;
2303 else gpar[paridx] += dpar/2;
2304 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2305 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2306
2307 Double_t h2 = 1./(2.*dpar);
2308 Double_t d0 = pintl4[0]-pintl1[0];
2309 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
2310 derivative[0] = h2*(4*d2 - d0)/3.;
2311 if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
2312
2313 d0 = pintl4[2]-pintl1[2];
2314 d2 = 2.*(pintl3[2]-pintl2[2]);
2315 derivative[2] = h2*(4*d2 - d0)/3.;
2316 if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
2317
2318 AliDebug(3,Form("\n+++ derivatives +++ \n"));
2319 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2320 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
2321
2322 return 0;
2323}
2324
2325//________________________________________________________________________________________________________
2326Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2327{
2328 /// Define local equation for current cluster in X and Z coor.
2329 /// and store them to memory
2330 /// return -1 in case of failure to build some equation
2331 /// 0 if no free global parameters were found but local eq is built
2332 /// 1 if both local and global eqs are built
2333 //
2334 // store first intersection point
2335 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
2336 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2337
2338 AliDebug(2,Form("Intersect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2339
2340 // calculate local derivatives numerically
2341 Bool_t zeroX = kTRUE;
2342 Bool_t zeroZ = kTRUE;
2343 //
2344 for (Int_t i=0; i<fNLocal; i++) {
2345 if (CalcDerivatives(i,kTRUE)) return -1;
2346 m.fDerLoc[i][kX] = fDerivativeLoc[i][0];
2347 m.fDerLoc[i][kZ] = fDerivativeLoc[i][2];
2348 if (zeroX) zeroX = IsZero(fDerivativeLoc[i][0]);
2349 if (zeroZ) zeroZ = IsZero(fDerivativeLoc[i][2]);
2350 }
2351 // for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
2352 //
2353 if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2354 if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
2355 //
2356 int ifill = 0;
2357 //
2358 AliITSAlignMille2Module* endModule = fCurrentModule;
2359 //
2360 zeroX = zeroZ = kTRUE;
2361 Bool_t dfDone[kNParCh];
2362 for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2363 m.fNModFilled = 0;
2364 //
2365 // special block for SDD derivatives
2366 Double_t jacobian[kNParChGeom];
2367 Int_t nmodTested = 0;
2368 //
2369 do {
2370 if (fCurrentModule->GetNParFree()==0) continue;
2371 nmodTested++;
2372 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2373 //
2374 if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2375 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2376 if (!dfDone[i]) {
2377 if (CalcDerivatives(i,kFALSE)) return -1;
2378 else {
2379 dfDone[i] = kTRUE;
2380 if (zeroX) zeroX = IsZero(fDerivativeGlo[i][0]);
2381 if (zeroZ) zeroZ = IsZero(fDerivativeGlo[i][2]);
2382 }
2383 }
2384 //
2385 m.fDerGlo[ifill][kX] = fDerivativeGlo[i][0];
2386 m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][2];
2387 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2388 }
2389 //
2390 // specific for special sensors
2391 if ( fCurrentModule->IsSDD() &&
2392 (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0 ||
2393 fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) ) {
2394 //
2395 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2396 // where V0 and T are the nominal drift velocity, time and time0
2397 // and the dT0 and dV are the corrections:
2398 // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2399 // dX/dV = dX/dxloc * dxloc/dV = dX/dxloc * (T-T0)
2400 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2401 //
2402 if (!dfDone[AliITSAlignMille2Module::kDOFT0] || !dfDone[AliITSAlignMille2Module::kDOFDV]) {
2403 //
2404 double dXdxlocsens=0., dZdxlocsens=0.;
2405 //
2406 // if the current module is the sensor itself and we work with local params, then
2407 // we can directly take dX/dxloc_sens dZ/dxloc_sens
2408 if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2409 if (!dfDone[AliITSAlignMille2Module::kDOFTX]) {
2410 CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2411 dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2412 }
2413 dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2414 dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2415 }
2416 //
2417 else { // need to perform some transformations
2418 // fetch the jacobian of the transformation from the sensors local frame to the frame
2419 // where the parameters are defined:
2420 // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2421 if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2422 AliITSAlignMille2Module::kDOFTX, jacobian);
2423 // Local: dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2424 else fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2425 AliITSAlignMille2Module::kDOFTX, jacobian);
2426 //
2427 for (int j=0;j<kNParChGeom;j++) {
2428 // need global derivative even if the j-th param is locked
2429 if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2430 dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2431 dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2432 }
2433 }
2434 //
2435 if (zeroX) zeroX = IsZero(dXdxlocsens);
2436 if (zeroZ) zeroZ = IsZero(dZdxlocsens);
2437 //
2438 double vdrift = GetVDriftSDD();
2439 double tdrift = GetTDriftSDD();
2440 //
2441 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2442 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2443 dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2444 //
2445 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0] = -dXdxlocsens*TMath::Sign(tdrift,vdrift);
2446 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2] = -dZdxlocsens*TMath::Sign(tdrift,vdrift);
2447 dfDone[AliITSAlignMille2Module::kDOFDV] = kTRUE;
2448 //
2449 }
2450 //
2451 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2452 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2453 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2454 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2455 }
2456 //
2457 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) {
2458 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0];
2459 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2];
2460 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV);
2461 }
2462 }
2463 //
2464 m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2465 } while( (fCurrentModule=fCurrentModule->GetParent()) );
2466 //
2467 if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2468 if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
2469 //
2470 // ok, can copy to m
2471 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2472 m.fMeas[kX] = fMeasLoc[0]-fPintLoc0[0];
2473 m.fSigma[kX] = fSigmaLoc[0];
2474 //
2475 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2476 m.fMeas[kZ] = fMeasLoc[2]-fPintLoc0[2];
2477 m.fSigma[kZ] = fSigmaLoc[2];
2478 //
2479 m.fNGlobFilled = ifill;
2480 fCurrentModule = endModule;
2481 //
2482 return Int_t(!zeroX && !zeroZ);
2483}
2484
2485//________________________________________________________________________________________________________
2486Int_t AliITSAlignMille2::AddLocalEquationTPA(Mille2Data &m)
2487{
2488 /// Define local equation for current cluster in X Y and Z coor.
2489 /// and store them to memory
2490 /// return -1 in case of failure to build some equation
2491 /// 0 if no free global parameters were found but local eq is built
2492 /// 1 if both local and global eqs are built
2493 //
2494 int curpoint = fCluster.GetUniqueID()-1;
2495 TGeoHMatrix *tempHMat = GetSensorCurrMatrixSID(fCurrentSensID);// fCurrentModule->GetSensitiveVolumeMatrix(fCluster.GetVolumeID());
2496 //
2497 fTPAFitter->GetDResDParams(&fDerivativeLoc[0][0], curpoint); // resid. derivatives over the track parameters
2498 for (Int_t i=fNLocal; i--;) tempHMat->MasterToLocalVect(fDerivativeLoc[i],m.fDerLoc[i]);
2499 //
2500 // derivatives over the global parameters ---------------------------------------->>>
2501 Double_t dRdP[3][3]; // derivative of local residuals vs local position
2502 Double_t dPdG[AliITSAlignMille2Module::kMaxParGeom][3]; // derivatives of local position vs global params
2503 fTPAFitter->GetDResDPos(&fDerivativeGlo[0][0], curpoint);
2504 for (int i=3;i--;) tempHMat->MasterToLocalVect(fDerivativeGlo[i],dRdP[i]);
2505 //
2506 UInt_t ifill=0, dfDone = 0;
2507 m.fNModFilled = 0;
2508 //
2509 AliITSAlignMille2Module* endModule = fCurrentModule;
2510 //
2511 do {
2512 if (fCurrentModule->GetNParFree()==0) continue;
2513 if (!fUseGlobalDelta) dfDone = 0; // for local deltas the derivatives at diff. levels are different
2514 Bool_t jacobOK = kFALSE;
2515 //
2516 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2517 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2518 //
2519 if (!TestWordBit(dfDone,i)) { // need to calculate new derivative
2520 if (!jacobOK) {fCurrentModule->CalcDerivDPosDPar(fCluster.GetVolumeID(),fMeasLoc,&dPdG[0][0]); jacobOK = kTRUE;}
2521 // dRes_j/dGlo_i = \sum_{k=1:3} dRes_j/dPos_k * dPos_k/dGlo_i
2522 fDerivativeGlo[i][kX] = dRdP[kX][kX]*dPdG[i][kX] + dRdP[kY][kX]*dPdG[i][kY] + dRdP[kZ][kX]*dPdG[i][kZ];
2523 fDerivativeGlo[i][kY] = dRdP[kX][kY]*dPdG[i][kX] + dRdP[kY][kY]*dPdG[i][kY] + dRdP[kZ][kY]*dPdG[i][kZ];
2524 fDerivativeGlo[i][kZ] = dRdP[kX][kZ]*dPdG[i][kX] + dRdP[kY][kZ]*dPdG[i][kY] + dRdP[kZ][kZ]*dPdG[i][kZ];
2525 SetWordBit(dfDone,i);
2526 }
2527 //
2528 m.fDerGlo[ifill][kX] = fDerivativeGlo[i][kX];
2529 m.fDerGlo[ifill][kY] = fDerivativeGlo[i][kY];
2530 m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][kZ];
2531 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2532 //
2533 }
2534 //
2535 if ( fCurrentModule->IsSDD() ) { // specific for SDD
2536 //
2537 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2538 // where V0 and T are the nominal drift velocity, time and time0
2539 // and the dT0 and dV are the corrections:
2540 // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2541 // dX/dV = dX/dxloc * dxloc/dV = dX/dxloc * (T-T0)
2542 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2543 //
2544 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2545 if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFT0)) {
2546 double vdrift = GetVDriftSDD();
2547 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX] = -dRdP[kX][kX]*vdrift;
2548 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY] = -dRdP[kX][kY]*vdrift;
2549 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ] = -dRdP[kX][kZ]*vdrift;
2550 SetWordBit(dfDone, AliITSAlignMille2Module::kDOFT0);
2551 }
2552 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX];
2553 m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY];
2554 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ];
2555 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2556 }
2557 //
2558 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) {
2559 if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFDV)) {
2560 double tdrift = TMath::Sign(GetTDriftSDD(), GetVDriftSDD());
2561 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kX] = dRdP[kX][kX]*tdrift;
2562 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kY] = dRdP[kX][kY]*tdrift;
2563 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kZ] = dRdP[kX][kZ]*tdrift;
2564 SetWordBit(dfDone, AliITSAlignMille2Module::kDOFDV);
2565 }
2566 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kX];
2567 m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kY];
2568 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][kZ];
2569 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV);
2570 }
2571 }
2572 //
2573 m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2574 } while( (fCurrentModule=fCurrentModule->GetParent()) );
2575 //
2576 // store first local residuals
2577 fTPAFitter->GetResiduals(fPintLoc , curpoint); // lab residuals
2578 for (int i=3;i--;) fPintLoc[i] = -fPintLoc[i];
2579 tempHMat->MasterToLocalVect(fPintLoc,m.fMeas); // local residuals
2580 m.fSigma[kX] = fSigmaLoc[kX];
2581 m.fSigma[kY] = fSigmaLoc[kY];
2582 m.fSigma[kZ] = fSigmaLoc[kZ];
2583 //
2584 m.fNGlobFilled = ifill;
2585 fCurrentModule = endModule;
2586 //
2587 return 1;
2588}
2589
2590//________________________________________________________________________________________________________
2591void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
2592{
2593 /// Set local equations with data stored in m
2594 /// return 0 if success
2595 //
2596 for (Int_t j=0; j<neq; j++) {
2597 //
2598 const Mille2Data &m = marr[j];
2599 //
2600 Bool_t filled = kFALSE;
2601 for (int ic=3;ic--;) {
2602 if (ic==kY && !fUseLocalYErr) continue;
2603 AliDebug(2,Form("setting local equation %c with fMeas=%.6f and fSigma=%.6f",fgkXYZ[ic],m.fMeas[ic], m.fSigma[ic]));
2604 Bool_t zero = kTRUE;
2605 for (int i=fNLocal; i--;) zero &= SetLocalDerivative( i, m.fDerLoc[i][ic] );
2606 for (int i=m.fNGlobFilled;i--;) zero &= SetGlobalDerivative( m.fParMilleID[i] , m.fDerGlo[i][ic] );
2607 if (zero) { AliInfo(Form("Skipping %c residual due to the zero derivatives!",fgkXYZ[ic])); continue; }
2608 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeas[ic], m.fSigma[ic]);
2609 filled = kTRUE;
2610 //
2611 }
2612 //
2613 if (filled) for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
2614 }
2615}
2616
2617//________________________________________________________________________________________________________
2618Int_t AliITSAlignMille2::GlobalFit()
2619{
2620 /// Call global fit; Global parameters are stored in parameters
2621 if (!fIsMilleInit) Init();
2622 //
2623 ApplyPreConstraints();
2624 int res = fMillepede->GlobalFit();
2625 AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
2626 if (res) {
2627 // fetch the parameters
2628 for (int imd=fNModules;imd--;) {
2629 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2630 int nprocp = 0;
2631 for (int ip=mod->GetNParTot();ip--;) {
2632 int idp = mod->GetParOffset(ip);
2633 if (idp<0) continue; // was not in the explicit fit
2634 mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
2635 mod->SetParErr(ip,fMillepede->GetFinalError(idp));
2636 int np = fMillepede->GetProcessedPoints(idp);
2637 if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
2638 }
2639 if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
2640 }
2641
2642 }
2643 ApplyPostConstraints();
2644 return res;
2645}
2646
2647//________________________________________________________________________________________________________
2648void AliITSAlignMille2::PrintGlobalParameters()
2649{
2650 /// Print global parameters
2651 if (!fIsMilleInit) {
2652 AliInfo("Millepede has not been initialized!");
2653 return;
2654 }
2655 fMillepede->PrintGlobalParameters();
2656}
2657
2658//________________________________________________________________________________________________________
2659Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
2660{
2661 // load definitions of supermodules from a root file
2662 // return 0 if success
2663
2664 TFile *smf=TFile::Open(sfile);
2665 if (!smf->IsOpen()) {
2666 AliInfo(Form("Cannot open supermodule file %s",sfile));
2667 return -1;
2668 }
2669
2670 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
2671 if (!sma) {
2672 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
2673 return -2;
2674 }
2675 Int_t nsma=sma->GetEntriesFast();
2676 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
2677 //
2678 // pepo200709
2679 Char_t st[2048];
2680 char symname[250];
2681 // end pepo200709
2682
2683 UShort_t volid;
2684 TGeoHMatrix m;
2685 //
2686 for (Int_t i=0; i<nsma; i++) {
2687 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
2688 volid=a->GetVolUID();
2689 strcpy(st,a->GetSymName());
2690 a->GetMatrix(m);
2691 //
2692 sscanf(st,"%s",symname);
2693 //
2694 // decode module list
2695 char *stp=strstr(st,"ModuleList:");
2696 if (!stp) return -3;
2697 stp += 11;
2698 int idx[2200];
2699 char spp[200]; int jp=0;
2700 char cl[20];
2701 strcpy(st,stp);
2702 int l=strlen(st);
2703 int j=0;
2704 int n=0;
2705 //
2706 while (j<=l) {
2707 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
2708 spp[jp]=0;
2709 jp=0;
2710 if (strlen(spp)) {
2711 int k=strcspn(spp,"-");
2712 if (k<int(strlen(spp))) { // c'e' il -
2713 strcpy(cl,&(spp[k+1]));
2714 spp[k]=0;
2715 int ifrom=atoi(spp); int ito=atoi(cl);
2716 for (int b=ifrom; b<=ito; b++) {
2717 idx[n]=b;
2718 n++;
2719 }
2720 }
2721 else { // numerillo singolo
2722 idx[n]=atoi(spp);
2723 n++;
2724 }
2725 }
2726 }
2727 else {
2728 spp[jp]=st[j];
2729 jp++;
2730 }
2731 j++;
2732 }
2733 UShort_t volidsv[2198];
2734 for (j=0;j<n;j++) {
2735 volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
2736 if (!volidsv[j]) {
2737 AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
2738 return -5;
2739 }
2740 }
2741 Int_t smindex=int(2198+volid-14336); // virtual index
2742 //
2743 fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
2744 //
2745 fNSuperModules++;
2746 }
2747 //
2748 smf->Close();
2749 //
2750 return 0;
2751}
2752
2753//________________________________________________________________________________________________________
2754void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
2755{
2756 // require that sum of modifications for the childs of this module is = val, i.e.
2757 // the internal corrections moves the module as a whole by fixed value (0 by default).
2758 // pattern is the bit pattern for the parameters to constrain
2759 //
2760 if (fIsMilleInit) {
2761 AliInfo("Millepede has been already initialized: no constrain may be added!");
2762 return;
2763 }
2764 if (!GetMilleModule(idm)->GetNChildren()) return;
2765 TString nm = "cstrSUMean";
2766 nm += GetNConstraints();
2767 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2768 idm,val,pattern);
2769 cstr->SetConstraintID(GetNConstraints());
2770 fConstraints.Add(cstr);
2771}
2772
2773//________________________________________________________________________________________________________
2774void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
2775{
2776 // require that median of the modifications for the childs of this module is = val, i.e.
2777 // the internal corrections moves the module as a whole by fixed value (0 by default)
2778 // module the outliers.
2779 // pattern is the bit pattern for the parameters to constrain
2780 // The difference between the mean and the median will be transfered to the parent
2781 if (fIsMilleInit) {
2782 AliInfo("Millepede has been already initialized: no constrain may be added!");
2783 return;
2784 }
2785 if (!GetMilleModule(idm)->GetNChildren()) return;
2786 TString nm = "cstrSUMed";
2787 nm += GetNConstraints();
2788 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2789 idm,val,pattern);
2790 cstr->SetConstraintID(GetNConstraints());
2791 fConstraints.Add(cstr);
2792}
2793
2794//________________________________________________________________________________________________________
2795void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
2796{
2797 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2798 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2799 // pattern is the bit pattern for the parameters to constrain
2800 //
2801 if (fIsMilleInit) {
2802 AliInfo("Millepede has been already initialized: no constrain may be added!");
2803 return;
2804 }
2805 TString nm = "cstrOMean";
2806 nm += GetNConstraints();
2807 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2808 -1,val,pattern);
2809 cstr->SetConstraintID(GetNConstraints());
2810 fConstraints.Add(cstr);
2811}
2812
2813//________________________________________________________________________________________________________
2814void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
2815{
2816 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2817 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2818 // pattern is the bit pattern for the parameters to constrain
2819 //
2820 if (fIsMilleInit) {
2821 AliInfo("Millepede has been already initialized: no constrain may be added!");
2822 return;
2823 }
2824 TString nm = "cstrOMed";
2825 nm += GetNConstraints();
2826 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2827 -1,val,pattern);
2828 cstr->SetConstraintID(GetNConstraints());
2829 fConstraints.Add(cstr);
2830}
2831
2832//________________________________________________________________________________________________________
2833void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
2834{
2835 // apply constraint on parameters in the local frame
2836 if (fIsMilleInit) {
2837 AliInfo("Millepede has been already initialized: no constrain may be added!");
2838 return;
2839 }
2840 AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
2841 cstr->SetConstraintID(GetNConstraints());
2842 fConstraints.Add(cstr);
2843}
2844
2845//________________________________________________________________________________________________________
2846void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
2847{
2848 // apply the constraint on the local corrections of a list of modules
2849 int nmod = cstr->GetNModules();
2850 double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
2851 //
2852 for (int imd=nmod;imd--;) {
2853 int modID = cstr->GetModuleID(imd);
2854 AliITSAlignMille2Module* mod = GetMilleModule(modID);
2855 ResetLocalEquation();
2856 int nadded = 0;
2857 double value = cstr->GetValue();
2858 double sigma = cstr->GetError();
2859 //
2860 // in case the reference (survey) deltas were imposed for Gaussian constraints
2861 // already accumulated corrections: they must be subtracted from the constraint value.
2862 if (IsConstraintWrtRef()) {
2863 //
2864 Double_t precal[AliITSAlignMille2Module::kMaxParTot];
2865 Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
2866 for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
2867 //
2868 // check if there was a reference delta provided for this module
2869 AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
2870 if (parref) parref->GetPars(refcal, refcal+3); // found reference delta
2871 //
2872 // extract already applied local corrections for this module
2873 if (fPrealignment) {
2874 //
2875 AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
2876 if (preo) {
2877 TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); // Delta_Glob * Delta_Glob_Par * M
2878 preo->GetMatrix(preMat); // Delta_Glob
2879 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
2880 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
2881 AliAlignObjParams algob;
2882 algob.SetMatrix(tmpMat);
2883 algob.GetPars(precal,precal+3); // local corrections for geometry
2884 }
2885 }
2886 //
2887 // subtract the contribution to constraint from precalibration
2888 for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
2889 //
2890 }
2891 //
2892 if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
2893 //
2894 for (int ipar=cstr->GetNCoeffs();ipar--;) {
2895 double coef = cstr->GetCoeff(ipar);
2896 if (IsZero(coef)) continue;
2897 //
2898 if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
2899 // we are working with local params or if the given param is not related to geometry,
2900 // apply the constraint directly
2901 int parPos = mod->GetParOffset(ipar);
2902 if (parPos<0) continue; // not in the fit
2903 fGlobalDerivatives[parPos] += coef;
2904 nadded++;
2905 }
2906 else { // we are working with global params, while the constraint is on local ones -> jacobian
2907 for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
2908 int parPos = mod->GetParOffset(jpar);
2909 if (parPos<0) continue;
2910 fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
2911 nadded++;
2912 }
2913 }
2914 }
2915 if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
2916 }
2917 //
2918}
2919
2920//________________________________________________________________________________________________________
2921void AliITSAlignMille2::ApplyPreConstraints()
2922{
2923 // apply constriants which cannot be imposed after the fit
2924 int nconstr = GetNConstraints();
2925 for (int i=0;i<nconstr;i++) {
2926 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2927 //
2928 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
2929 ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
2930 continue;
2931 }
2932 //
2933 if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
2934 //
2935 if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
2936 // apply constraint on the mean's before the fit
2937 int imd = cstr->GetModuleID();
2938 if (imd>=0) {
2939 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2940 UInt_t pattern = 0;
2941 for (int ipar=mod->GetNParTot();ipar--;) {
2942 if (!cstr->IncludesParam(ipar)) continue;
2943 if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
2944 pattern |= 0x1<<ipar;
2945 cstr->SetApplied(ipar);
2946 }
2947 ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
2948 //
2949 }
2950 else if (!PseudoParentsAllowed()) {
2951 ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
2952 cstr->SetApplied(-1);
2953 }
2954 }
2955}
2956
2957//________________________________________________________________________________________________________
2958void AliITSAlignMille2::ApplyPostConstraints()
2959{
2960 // apply constraints which can be imposed after the fit
2961 int nconstr = GetNConstraints();
2962 Bool_t convGlo = kFALSE;
2963 // check if there is something to do
2964 int ntodo = 0;
2965 for (int i=0;i<nconstr;i++) {
2966 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2967 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2968 if (cstr->GetRemainingPattern() == 0) continue;
2969 ntodo++;
2970 }
2971 if (!ntodo) return;
2972 //
2973 if (!fUseGlobalDelta) { // need to convert to global params
2974 ConvertParamsToGlobal();
2975 convGlo = kTRUE;
2976 }
2977 //
2978 for (int i=0;i<nconstr;i++) {
2979 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2980 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2981 //
2982 int imd = cstr->GetModuleID();
2983 //
2984 if (imd>=0) {
2985 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2986 UInt_t pattern = 0;
2987 for (int ipar=mod->GetNParTot();ipar--;) {
2988 if (cstr->IsApplied(ipar)) continue;
2989 if (!cstr->IncludesParam(ipar)) continue;
2990 if (!mod->IsFreeDOF(ipar)) continue; // parameter is fixed, will not apply constraint
2991 pattern |= 0x1<<ipar;
2992 cstr->SetApplied(ipar);
2993 }
2994 if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
2995 //
2996 }
2997 else if (PseudoParentsAllowed()) {
2998 UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
2999 PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
3000 cstr->SetApplied(-1);
3001 }
3002 }
3003 // if there was a conversion, rewind it
3004 if (convGlo) ConvertParamsToLocal();
3005 //
3006}
3007
3008//________________________________________________________________________________________________________
3009void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
3010{
3011 // require that sum of modifications for the childs of this module is = val, i.e.
3012 // the internal corrections moves the module as a whole by fixed value (0 by default).
3013 // pattern is the bit pattern for the parameters to constrain
3014 //
3015 //
3016 AliITSAlignMille2Module* mod = GetMilleModule(idm);
3017 //
3018 for (int ip=0;ip<kNParCh;ip++) {
3019 if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
3020 ResetLocalEquation();
3021 int nadd = 0;
3022 for (int ich=mod->GetNChildren();ich--;) {
3023 int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
3024 if (idpar<0) continue;
3025 fGlobalDerivatives[idpar] = 1.0;
3026 nadd++;
3027 }
3028 //
3029 if (nadd>0) {
3030 AddConstraint(fGlobalDerivatives,val);
3031 AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
3032 }
3033 }
3034 //
3035}
3036
3037//________________________________________________________________________________________________________
3038void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
3039{
3040 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3041 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3042 // pattern is the bit pattern for the parameters to constrain
3043 //
3044 for (int ip=0;ip<kNParCh;ip++) {
3045 //
3046 if ( !((pattern>>ip)&0x1) ) continue;
3047 ResetLocalEquation();
3048 int nadd = 0;
3049 for (int imd=fNModules;imd--;) {
3050 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3051 if (mod->GetParent()) continue; // this is not an orphan
3052 int idpar = mod->GetParOffset(ip);
3053 if (idpar<0) continue;
3054 fGlobalDerivatives[idpar] = 1.0;
3055 nadd++;
3056 }
3057 if (nadd>0) {
3058 AddConstraint(fGlobalDerivatives,val);
3059 AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
3060 }
3061 }
3062 //
3063 //
3064}
3065
3066//________________________________________________________________________________________________________
3067void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
3068{
3069 // require that median or mean of the modifications for the childs of this module is = val, i.e.
3070 // the internal corrections moves the module as a whole by fixed value (0 by default)
3071 // module the outliers.
3072 // pattern is the bit pattern for the parameters to constrain
3073 // The difference between the mean and the median will be transfered to the parent
3074 //
3075 AliITSAlignMille2Module* parent = GetMilleModule(idm);
3076 int nc = parent->GetNChildren();
3077 //
3078 double *tmpArr = new double[nc];
3079 //
3080 for (int ip=0;ip<kNParCh;ip++) {
3081 int npc = 0;
3082 if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
3083 // compute the mean and median of the deltas
3084 int nfree = 0;
3085 for (int ich=nc;ich--;) {
3086 AliITSAlignMille2Module* child = parent->GetChild(ich);
3087 // if (!child->IsFreeDOF(ip)) continue;
3088 tmpArr[nfree++] = child->GetParVal(ip);
3089 }
3090 double median=0,mean=0;
3091 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3092 mean += tmpArr[ic0];
3093 for (int ic1=ic0+1;ic1<nfree;ic1++)
3094 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3095 }
3096 //
3097 int kmed = nfree/2;
3098 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3099 if (nfree>0) mean /= nfree;
3100 //
3101 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3102 //
3103 for (int ich=nc;ich--;) {
3104 AliITSAlignMille2Module* child = parent->GetChild(ich);
3105 // if (!child->IsFreeDOF(ip)) continue;
3106 child->SetParVal(ip, child->GetParVal(ip) + shift);
3107 npc++;
3108 }
3109 //
3110 parent->SetParVal(ip, parent->GetParVal(ip) - shift);
3111 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d children of module %d: %s",
3112 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3113 ip,npc,idm,parent->GetName()));
3114 }
3115 delete[] tmpArr;
3116 //
3117 //
3118}
3119
3120//________________________________________________________________________________________________________
3121void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
3122{
3123 // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
3124 // the corrections moves the whole setup by fixed value (0 by default).
3125 // pattern is the bit pattern for the parameters to constrain
3126 //
3127 int nc = fNModules;
3128 //
3129 int norph = 0;
3130 for (int ich=nc;ich--;) if (!GetMilleModule(ich)->GetParent()) norph ++;
3131 if (!norph) return;
3132 double *tmpArr = new double[norph];
3133 //
3134 for (int ip=0;ip<kNParCh;ip++) {
3135 int npc = 0;
3136 if ( !((pattern>>ip)&0x1)) continue;
3137 // compute the mean and median of the deltas
3138 int nfree = 0;
3139 for (int ich=nc;ich--;) {
3140 AliITSAlignMille2Module* child = GetMilleModule(ich);
3141 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3142 if (child->GetParent()) continue;
3143 tmpArr[nfree++] = child->GetParVal(ip);
3144 }
3145 double median=0,mean=0;
3146 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3147 mean += tmpArr[ic0];
3148 for (int ic1=ic0+1;ic1<nfree;ic1++)
3149 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3150 }
3151 //
3152 int kmed = nfree/2;
3153 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3154 if (nfree>0) mean /= nfree;
3155 //
3156 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3157 //
3158 for (int ich=nc;ich--;) {
3159 AliITSAlignMille2Module* child = GetMilleModule(ich);
3160 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
3161 if (child->GetParent()) continue;
3162 child->SetParVal(ip, child->GetParVal(ip) + shift);
3163 npc++;
3164 }
3165 //
3166 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d orphan modules",
3167 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3168 ip,npc));
3169 }
3170 delete[] tmpArr;
3171 //
3172}
3173
3174//________________________________________________________________________________________________________
3175Bool_t AliITSAlignMille2::IsParModConstrained(const AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
3176{
3177 // check if par of the module participates in some constraint, and set the flag for their types
3178 meanmed = gaussian = kFALSE;
3179 //
3180 if ( mod->IsParConstrained(par) ) gaussian = kTRUE; // direct constraint on this param
3181 //
3182 for (int icstr=GetNConstraints();icstr--;) {
3183 AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
3184 //
3185 if (!cstr->IncludesModPar(mod,par)) continue;
3186 if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
3187 else meanmed = kTRUE;
3188 //
3189 if (meanmed && gaussian) break; // no sense to check further
3190 }
3191 //
3192 return meanmed||gaussian;
3193}
3194
3195//________________________________________________________________________________________________________
3196Bool_t AliITSAlignMille2::IsParModFamilyVaried(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3197{
3198 // check if parameter par is varied for this module or its children up to the level depth
3199 if (depth<0) return kFALSE;
3200 if (mod->GetParOffset(par)>=0) return kTRUE;
3201 for (int icld=mod->GetNChildren();icld--;) {
3202 AliITSAlignMille2Module* child = mod->GetChild(icld);
3203 if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
3204 }
3205 return kFALSE;
3206 //
3207}
3208
3209/*
3210//________________________________________________________________________________________________________
3211Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3212{
3213 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3214 if (depth<0) return kTRUE;
3215 for (int icld=mod->GetNChildren();icld--;) {
3216 AliITSAlignMille2Module* child = mod->GetChild(icld);
3217 //if (child->GetParOffset(par)<0) continue; // fixed
3218 Bool_t cstMM=kFALSE,cstGS=kFALSE;
3219 // does this child have gaussian constraint ?
3220 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3221 // check its children
3222 if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
3223 }
3224 return kFALSE;
3225 //
3226}
3227*/
3228
3229//________________________________________________________________________________________________________
3230Bool_t AliITSAlignMille2::IsParFamilyFree(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3231{
3232 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3233 if (depth<0) return kFALSE;
3234 for (int icld=mod->GetNChildren();icld--;) {
3235 AliITSAlignMille2Module* child = mod->GetChild(icld);
3236 //if (child->GetParOffset(par)<0) continue; // fixed
3237 Bool_t cstMM=kFALSE,cstGS=kFALSE;
3238 // does this child have gaussian constraint ?
3239 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3240 // check its children
3241 if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
3242 }
3243 return kFALSE;
3244 //
3245}
3246
3247//________________________________________________________________________________________________________
3248Double_t AliITSAlignMille2::GetTDriftSDD() const
3249{
3250 // obtain drift time corrected for t0
3251 double t = fCluster.GetDriftTime();
3252 return t - fDriftTime0[ fCluster.GetUniqueID()-1 ];
3253}
3254
3255//________________________________________________________________________________________________________
3256Double_t AliITSAlignMille2::GetVDriftSDD() const
3257{
3258 // obtain corrected drift speed
3259 return fDriftSpeed[ fCluster.GetUniqueID()-1 ];
3260}
3261
3262//________________________________________________________________________________________________________
3263Bool_t AliITSAlignMille2::FixedOrphans() const
3264{
3265 // are there fixed modules with no parent (normally in such a case
3266 // the constraints on the orphans should not be applied
3267 if (!IsConfigured()) {
3268 AliInfo("Still not configured");
3269 return kFALSE;
3270 }
3271 for (int i=0;i<fNModules;i++) {
3272 AliITSAlignMille2Module* md = GetMilleModule(i);
3273 if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
3274 }
3275 return kFALSE;
3276}
3277
3278//________________________________________________________________________________________________________
3279void AliITSAlignMille2::ConvertParamsToGlobal()
3280{
3281 // convert params in local frame to global one
3282 double pars[AliITSAlignMille2Module::kMaxParGeom];
3283 for (int imd=fNModules;imd--;) {
3284 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3285 if (mod->GeomParamsGlobal()) continue;
3286 mod->GetGeomParamsGlo(pars);
3287 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3288 mod->SetGeomParamsGlobal(kTRUE);
3289 }
3290}
3291
3292//________________________________________________________________________________________________________
3293void AliITSAlignMille2::ConvertParamsToLocal()
3294{
3295 // convert params in global frame to local one
3296 double pars[AliITSAlignMille2Module::kMaxParGeom];
3297 for (int imd=fNModules;imd--;) {
3298 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3299 if (!mod->GeomParamsGlobal()) continue;
3300 mod->GetGeomParamsLoc(pars);
3301 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3302 mod->SetGeomParamsGlobal(kFALSE);
3303 }
3304}
3305
3306//________________________________________________________________________________________________________
3307void AliITSAlignMille2::SetBField(Double_t b)
3308{
3309 // set Bz value
3310 if (IsZero(b,1e-5)) {
3311 fBField = 0.0;
3312 fBOn = kFALSE;
3313 fNLocal = 4;
3314 }
3315 else {
3316 fBField = b;
3317 fBOn = kTRUE;
3318 fNLocal = 5; // helices
3319 }
3320}
3321
3322//________________________________________________________________________________________________________
3323Int_t AliITSAlignMille2::ProcessUserInfo(TList* userInfo)
3324{
3325 // extract calibration information used for TrackPointArray creation from run info
3326 //
3327 if (!userInfo) { AliInfo("No UserInfo is provided"); return 0;}
3328 //
3329 TMap *cdbMap=0;
3330 TList* cdbList=0;
3331 TObjString *objStr,*keyStr;
3332 AliCDBManager* man = AliCDBManager::Instance();
3333 //
3334 int run = userInfo->GetUniqueID();
3335 AliInfo(Form("UserInfo corresponds to run#%d",run));
3336 cdbMap = (TMap*)userInfo->FindObject("cdbMap");
3337 if (!cdbMap) {AliInfo("No CDB Map found in UserInfo");}
3338 else {
3339 if ((objStr=(TObjString*)cdbMap->GetValue("default"))) { // first set default CDB path
3340 fDefCDBpath = objStr->GetString();
3341 if (fDefCDBpath.BeginsWith("raw://")) fDefCDBpath = "raw://";
3342 AliInfo(Form("Default CDB Storage from UserInfo: %s",fDefCDBpath.Data()));
3343 }
3344 man->SetDefaultStorage( fDefCDBpath.Data() ); // this may be overriden later by configuration file
3345 man->SetRun(run);
3346 //
3347 // set specific paths relevant for alignment
3348 TIter itMap(cdbMap);
3349 while( (keyStr=(TObjString*)itMap.Next()) ) {
3350 TString keyS = keyStr->GetString();
3351 if ( keyS == "default" ) continue;
3352 man->SetSpecificStorage( keyS.Data(), cdbMap->GetValue(keyS)->GetName() );
3353 }
3354 }
3355 //
3356 cdbList = (TList*)userInfo->FindObject("cdbList");
3357 if (!cdbList) {AliInfo("No CDB List found in UserInfo");}
3358 else {
3359 // Deltas used for TrackPointArray production
3360 TIter itList(cdbList);
3361 while( (objStr=(TObjString*)itList.Next()) )
3362 if (objStr->GetString().Contains("ITS/Align/Data")) {
3363 fInitDeltaPath = objStr->GetString();
3364 AliInfo(Form("Production Misalignment from UserInfo: %s",fInitDeltaPath.Data()));
3365 break;
3366 }
3367 // SDD response (time0 and drift speed correction) used for TrackPointArray production
3368 itList.Reset();
3369 while( (objStr=(TObjString*)itList.Next()) )
3370 if (objStr->GetString().Contains("ITS/Calib/RespSDD")) {
3371 fInitSDDRespPath = objStr->GetString();
3372 AliInfo(Form("Production SDD Response from UserInfo: %s",fInitSDDRespPath.Data()));
3373 break;
3374 }
3375 //
3376 }
3377 //
3378 objStr = (TObjString*)userInfo->FindObject("BzkGauss");
3379 if (objStr) {
3380 SetBField( objStr->GetString().Atof() );
3381 AliInfo(Form("Magentic field from UserInfo: %+.2e",GetBField()));
3382 }
3383 return 0;
3384}
3385
3386//________________________________________________________________________________________________________
3387Int_t AliITSAlignMille2::LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
3388{
3389 if (path.IsNull()) return 0;
3390 //
abd7ef79 3391 AliCDBEntry *entry = 0;
6be22b3f 3392 resp = 0;
3393 while(1) {
3394 if (path.BeginsWith("path: ")) { // must load from OCDB
3395 AliCDBId* cdbId = AliCDBId::MakeFromString( path.Data() );
abd7ef79 3396 entry = AliCDBManager::Instance()->Get( *cdbId );
6be22b3f 3397 delete cdbId;
3398 if (!entry) break;
3399 resp = (AliITSresponseSDD*) entry->GetObject();
3400 entry->SetObject(NULL);
3401 entry->SetOwner(kTRUE);
3402 delete entry;
3403 break;
3404 }
3405 //
3406 if (gSystem->AccessPathName(path.Data())) break;
3407 TFile* precf = TFile::Open(path.Data());
abd7ef79 3408 if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
3409 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
3410 resp = (AliITSresponseSDD*) entry->GetObject();
3411 if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
3412 else resp = 0;
3413 entry->SetOwner(kTRUE);
3414 delete entry;
3415 }
3416 //
6be22b3f 3417 precf->Close();
3418 delete precf;
3419 break;
3420 }
3421 //
3422 if (!resp) {AliError(Form("Failed to load SDD response from %s",path.Data())); return -1;}
3423 return 0;
3424}
3425
3426//________________________________________________________________________________________________________
3427Int_t AliITSAlignMille2::LoadDeltas(TString& path, TClonesArray *&arr)
3428{
3429 if (path.IsNull()) return 0;
3430 //
abd7ef79 3431 AliCDBEntry *entry = 0;
6be22b3f 3432 arr = 0;
3433 while(1) {
3434 if (path.BeginsWith("path: ")) { // must load from OCDB
3435 AliCDBId *cdbId = AliCDBId::MakeFromString( path.Data() );
abd7ef79 3436 entry = AliCDBManager::Instance()->Get( *cdbId );
6be22b3f 3437 delete cdbId;
3438 if (!entry) break;
3439 arr = (TClonesArray*) entry->GetObject();
3440 entry->SetObject(NULL);
3441 entry->SetOwner(kTRUE);
3442 delete entry;
3443 break;
3444 }
3445 //
3446 if (gSystem->AccessPathName(path.Data())) break;
3447 TFile* precf = TFile::Open(path.Data());
abd7ef79 3448 if (precf->FindKey("ITSAlignObjs")) arr = (TClonesArray*)precf->Get("ITSAlignObjs");
3449 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
3450 arr = (TClonesArray*) entry->GetObject();
3451 if (arr && arr->InheritsFrom(TClonesArray::Class())) entry->SetObject(NULL);
3452 else arr = 0;
3453 entry->SetOwner(kTRUE);
3454 delete entry;
3455 }
6be22b3f 3456 precf->Close();
3457 delete precf;
3458 break;
3459 }
3460 //
3461 if (!arr) {AliError(Form("Failed to load Deltas from %s",path.Data())); return -1;}
3462 return 0;
3463}
3464
3465//________________________________________________________________________________________________________
abd7ef79 3466Int_t AliITSAlignMille2::CacheMatricesCurr()
6be22b3f 3467{
3468 // build arrays for the fast access to sensor matrices from their sensor ID
3469 //
3470 TGeoHMatrix mdel;
abd7ef79 3471 AliInfo("Building sensors current matrices cache");
6be22b3f 3472 //
6be22b3f 3473 fCacheMatrixCurr.Delete();
6be22b3f 3474 for (int idx=0;idx<=kMaxITSSensID;idx++) {
3475 int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
6be22b3f 3476 TGeoHMatrix *mcurr = new TGeoHMatrix();
abd7ef79 3477 AliITSAlignMille2Module::SensVolMatrix(volID, mcurr);
6be22b3f 3478 fCacheMatrixCurr.AddAtAndExpand(mcurr,idx);
3479 //
3480 }
3481 //
6be22b3f 3482 fCacheMatrixCurr.SetOwner(kTRUE);
abd7ef79 3483 return 0;
3484}
3485
3486//________________________________________________________________________________________________________
3487Int_t AliITSAlignMille2::CacheMatricesOrig()
3488{
3489 // build arrays for the fast access to sensor original matrices (used for production)
3490 //
3491 TGeoHMatrix mdel;
3492 AliInfo("Building sensors original matrices cache");
3493 //
3494 fCacheMatrixOrig.Delete();
3495 if (!fInitDeltaPath.IsNull()) {
3496 if (LoadDeltas(fInitDeltaPath,fPrealignment) || ApplyToGeometry())
3497 { AliError("Failed to load/apply initial deltas used to produce points"); return -1;}
3498 }
3499 //
3500 for (int idx=0;idx<=kMaxITSSensID;idx++) {
3501 int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
3502 TGeoHMatrix *morig = new TGeoHMatrix();
3503 if (fUsePreAlignment) AliITSAlignMille2Module::SensVolMatrix(volID,morig);
3504 else AliITSAlignMille2Module::SensVolOrigGlobalMatrix(volID,morig);
3505 fCacheMatrixOrig.AddAtAndExpand(morig,idx);
3506 }
3507 //
3508 fCacheMatrixOrig.SetOwner(kTRUE);
3509 if (fUsePreAlignment) { // the initial deltas were temporary attached to prealignment array, clean and reinitialize geometry
3510 delete fPrealignment;
3511 fPrealignment = 0;
3512 fUsePreAlignment = 0;
3513 InitGeometry();
3514 }
6be22b3f 3515 return 0;
3516}
3517
3518//________________________________________________________________________________________________________
3519void AliITSAlignMille2::ConstrainHelixFitPT(Int_t q,Double_t pt,Double_t pterr)
3520{
3521 // constrain q and pT of the helical fit of the track (should be set before process.track)
3522 //
3523 fConstrCharge = q==0 ? q:TMath::Sign(1,q);
3524 fConstrPT = pt;
3525 fConstrPTErr = pterr;
3526}
3527
3528//________________________________________________________________________________________________________
3529void AliITSAlignMille2::ConstrainHelixFitCurv(Int_t q,Double_t crv,Double_t crverr)
3530{
3531 // constrain charge and curvature of the helical fit of the track (should be set before process.track)
3532 //
3533 const double kCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
3534
3535 fConstrCharge = q==0 ? q:TMath::Sign(1,q);
3536 if (crv<0 || IsZero(crv)) {
3537 fConstrPT = -1;
3538 fConstrPTErr = -1;
3539 }
3540 else {
3541 fConstrPT = 1./crv*fBField*kCQConv;
3542 fConstrPTErr =fConstrPT/crv*crverr;
3543 }
3544}
3545
3546//________________________________________________________________________________________________________
3547TClonesArray* AliITSAlignMille2::CreateDeltas()
3548{
3549 // Create \Deltas for every explicitly or implicitly (via non-alignable volumes) varied
3550 // or prealigned module.
3551 // If the module has inded J in the hierarchy of alignable volumes (0 - the top, most
3552 // coarse level), then its Delta is expressed via MP2 \deltas (in global frame) and
3553 // prealignment \DeltaP's as:
3554 // \Delta_J = Y X Y^-1
3555 // where X = \delta_J * \DeltaP_J
3556 // Y = Prod_{K=0,J-1} \delta_K
3557 // Note that \delta_L accounts not only for its own correction but also of all non-alignable
3558 // modules in the hierarchy chain from L up to the closest alignable:
3559 // while (parent && !parent->IsAlignable()) {
3560 // \delta_L->MultiplyLeft( \delta_parent );
3561 // parent = parent->GetParent();
3562 // }
3563 //
3564 Bool_t convLoc = kFALSE;
3565 if (!GetUseGlobalDelta()) {
3566 ConvertParamsToGlobal();
3567 convLoc = kTRUE;
3568 }
3569 //
3570 AliAlignObjParams tempAlignObj;
3571 TGeoHMatrix tempMatX,tempMatY,tempMat1;
3572 //
3573 TClonesArray *array = new TClonesArray("AliAlignObjParams",10);
3574 TClonesArray &alobj = *array;
3575 int idx = 0;
3576 //
3577 TGeoManager* geoManager = AliGeomManager::GetGeometry();
3578 int nalgtot = geoManager->GetNAlignable();
3579 //
3580 for (int ialg=0;ialg<nalgtot;ialg++) { // loop over all alignable entries
3581 //
3582 const char* algname = geoManager->GetAlignableEntry(ialg)->GetName();
3583 //
3584 AliITSAlignMille2Module* md = GetMilleModuleBySymName(algname); // explicitly varied?
3585 AliITSAlignMille2Module* parent = md ? md->GetParent(): GetMilleModuleIfContained(algname);
3586 AliAlignObjParams* preob = GetPrealignedObject(algname); // was it prealigned ?
3587 //
3588 if (!preob && !md && (!parent || parent->IsAlignable())) continue; // noting to do
3589 //
3590 // create matrix X (see comment) ------------------------------------------------->>>
3591 // start from unity matrix
3592 tempMatX.Clear();
3593 if (preob) { // account prealigngment
3594 preob->GetMatrix(tempMat1);
3595 tempMatX.MultiplyLeft(&tempMat1);
3596 }
3597 //
3598 if (md) {
3599 tempAlignObj.SetTranslation( md->GetParVal(0),md->GetParVal(1),md->GetParVal(2));
3600 tempAlignObj.SetRotation( md->GetParVal(3),md->GetParVal(4),md->GetParVal(5));
3601 tempAlignObj.GetMatrix(tempMat1);
3602 tempMatX.MultiplyLeft(&tempMat1); // acount correction to varied module
3603 }
3604 //
3605 // the corrections to all non-alignable modules from current on
3606 // till first alignable should add up to its matrix
3607 while (parent && !parent->IsAlignable()) {
3608 tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
3609 tempAlignObj.SetRotation( parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
3610 tempAlignObj.GetMatrix(tempMat1);
3611 tempMatX.MultiplyLeft(&tempMat1); // add matrix of non-alignable module
3612 parent = parent->GetParent();
3613 }
3614 // create matrix X (see comment) ------------------------------------------------<<<
3615 //
3616 // create matrix Y (see comment) ------------------------------------------------>>>
3617 // start from unity matrix
3618 tempMatY.Clear();
3619 while ( parent ) {
3620 tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
3621 tempAlignObj.SetRotation( parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
3622 tempAlignObj.GetMatrix(tempMat1);
3623 tempMatY.MultiplyLeft(&tempMat1);
3624 parent = parent->GetParent();
3625 }
3626 // create matrix Y (see comment) ------------------------------------------------<<<
3627 //
3628 tempMatX.MultiplyLeft(&tempMatY);
3629 tempMatX.Multiply(&tempMatY.Inverse());
3630 //
3631 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname);
3632 new(alobj[idx++]) AliAlignObjParams(algname,vid,tempMatX,kTRUE);
3633 //
3634 }
3635 //
3636 if (convLoc) ConvertParamsToLocal();
3637 //
3638 return array;
3639 //
3640}
3641
3642//_______________________________________________________________________________________
3643AliITSresponseSDD* AliITSAlignMille2::CreateSDDResponse()
3644{
3645 // create object with SDD repsonse (t0 and vdrift corrections) accounting for
3646 // eventual precalibration
3647 //
3648 // if there was a precalibration provided, copy it to new arrray
3649 AliITSresponseSDD *precal = GetSDDPrecalibration();
3650 if (!precal) precal = GetSDDInit();
3651 AliITSresponseSDD *calibSDD = new AliITSresponseSDD();
3652 //
3653 for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) {
3654 calibSDD->SetModuleTimeZero(ind, precal? precal->GetTimeZero(ind) : 0.);
3655 calibSDD->SetDeltaVDrift(ind, precal? precal->GetDeltaVDrift(ind) : 0.);
3656 }
3657 //
3658 Bool_t save = kFALSE;
3659 for (int imd=GetNModules();imd--;) {
3660 AliITSAlignMille2Module* md = GetMilleModule(imd);
3661 if (!md->IsSDD()) continue;
3662 if (md->IsFreeDOF(AliITSAlignMille2Module::kDOFT0) ||
3663 md->IsFreeDOF(AliITSAlignMille2Module::kDOFDV)) save = kTRUE;
3664 //
3665 for (int is=0;is<md->GetNSensitiveVolumes();is++) {
3666 int ind = md->GetSensVolIndex(is);
3667 float t0 = calibSDD->GetTimeZero(ind) + md->GetParVal(AliITSAlignMille2Module::kDOFT0);
3668 float dv = calibSDD->GetDeltaVDrift(ind) + md->GetParVal(AliITSAlignMille2Module::kDOFDV);
3669 //
3670 calibSDD->SetModuleTimeZero(ind, t0);
3671 calibSDD->SetDeltaVDrift(ind, dv);
3672 }
3673 }
3674 //
3675 if (!save) {
3676 AliInfo("No free parameters for SDD calibration, nothing to save");
3677 delete calibSDD;
3678 calibSDD = 0;
3679 }
3680 //
3681 return calibSDD;
3682}