]> 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
Read data format and sampling frequency from CDH in SDD DAs (F. Prino)
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille2.cxx
CommitLineData
7b85e477 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
057304f3 16/* $Id$ */
7b85e477 17
18//-----------------------------------------------------------------------------
6526a72c 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
7b85e477 28//-----------------------------------------------------------------------------
29
7b85e477 30#include <TFile.h>
31#include <TClonesArray.h>
7b85e477 32#include <TMath.h>
7b85e477 33#include <TVirtualFitter.h>
6526a72c 34#include <TGeoManager.h>
45993b23 35#include <TArrayI.h>
36#include <TSystem.h>
7b85e477 37#include "AliITSAlignMille2.h"
38#include "AliITSgeomTGeo.h"
39#include "AliGeomManager.h"
40#include "AliMillePede2.h"
41#include "AliTrackPointArray.h"
42#include "AliAlignObjParams.h"
43#include "AliLog.h"
7b85e477 44#include "AliTrackFitterRieman.h"
45993b23 45#include "AliITSAlignMille2Constraint.h"
46#include "AliITSAlignMille2ConstrArray.h"
47#include "AliITSresponseSDD.h"
6526a72c 48
7b85e477 49ClassImp(AliITSAlignMille2)
6526a72c 50
18986853 51const Char_t* AliITSAlignMille2::kRecKeys[] = {
52 "GEOMETRY_FILE",
53 "SUPERMODULE_FILE",
54 "CONSTRAINTS_REFERENCE_FILE",
55 "PREALIGNMENT_FILE",
56 "PRECALIBSDD_FILE",
57 "INITCALBSDD_FILE",
58 "SET_GLOBAL_DELTAS",
59 "CONSTRAINT_LOCAL",
60 "MODULE_VOLUID",
61 "MODULE_INDEX",
62 "SET_PSEUDO_PARENTS",
63 "SET_TRACK_FIT_METHOD",
64 "SET_MINPNT_TRA",
65 "SET_NSTDDEV",
66 "SET_RESCUT_INIT",
67 "SET_RESCUT_OTHER",
68 "SET_LOCALSIGMAFACTOR",
69 "SET_STARTFAC",
70 "SET_B_FIELD",
71 "SET_SPARSE_MATRIX",
72 "REQUIRE_POINT",
73 "CONSTRAINT_ORPHANS",
74 "CONSTRAINT_SUBUNITS",
75 "APPLY_CONSTRAINT"
76};
77
6526a72c 78
79//========================================================================================================
7b85e477 80
81AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;
82Int_t AliITSAlignMille2::fgInstanceID = 0;
83
6526a72c 84//________________________________________________________________________________________________________
85AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename )
7b85e477 86: TObject(),
87 fMillepede(0),
88 fStartFac(16.),
89 fResCutInitial(100.),
90 fResCut(100.),
91 fNGlobal(0),
92 fNLocal(4),
93 fNStdDev(3),
94 fIsMilleInit(kFALSE),
6526a72c 95 fAllowPseudoParents(kFALSE),
96 //
7b85e477 97 fCurrentModule(0),
98 fTrack(0),
6526a72c 99 fTrackBuff(0),
7b85e477 100 fCluster(),
101 fGlobalDerivatives(0),
6526a72c 102 //
7b85e477 103 fMinNPtsPerTrack(3),
104 fInitTrackParamsMeth(1),
105 fTotBadLocEqPoints(0),
106 fRieman(0),
107 //
6526a72c 108 fConstraints(0),
109 //
7b85e477 110 fUseGlobalDelta(kFALSE),
111 fRequirePoints(kFALSE),
112 fTempExcludedModule(-1),
113 //
114 fGeometryFileName("geometry.root"),
115 fPreAlignmentFileName(""),
6526a72c 116 fConstrRefFileName(""),
7b85e477 117 fGeoManager(0),
118 fIsConfigured(kFALSE),
119 fPreAlignQF(0),
120//
6526a72c 121 fCorrectSDD(0),
122 fInitialRecSDD(0),
7b85e477 123 fPrealignment(0),
6526a72c 124 fConstrRef(0),
7b85e477 125 fMilleModule(2),
126 fSuperModule(2),
127 fNModules(0),
128 fNSuperModules(0),
129 fUsePreAlignment(kFALSE),
7b85e477 130 fBOn(kFALSE),
131 fBField(0.0),
6526a72c 132 fBug(0),
133 fMilleVersion(2)
7b85e477 134{
135 /// main constructor that takes input from configuration file
7b85e477 136 for (int i=3;i--;) fSigmaFactor[i] = 1.0;
137 //
138 // new RS
139 for (Int_t i=0; i<6; i++) {
140 fNReqLayUp[i]=0;
141 fNReqLayDown[i]=0;
142 fNReqLay[i]=0;
143 }
144 for (Int_t i=0; i<3; i++) {
145 fNReqDetUp[i]=0;
146 fNReqDetDown[i]=0;
147 fNReqDet[i]=0;
148 }
149 //
150 Int_t lc=LoadConfig(configFilename);
151 if (lc) {
6526a72c 152 AliError(Form("Error %d loading configuration from %s",lc,configFilename));
153 exit(1);
7b85e477 154 }
155 //
6526a72c 156 fMillepede = new AliMillePede2();
7b85e477 157 fgInstance = this;
158 fgInstanceID++;
159 //
160}
161
6526a72c 162//________________________________________________________________________________________________________
163AliITSAlignMille2::~AliITSAlignMille2()
164{
7b85e477 165 /// Destructor
6526a72c 166 if (fMillepede) delete fMillepede; fMillepede = 0;
167 if (fGlobalDerivatives) delete[] fGlobalDerivatives; fGlobalDerivatives = 0;
168 if (fRieman) delete fRieman; fRieman = 0;
169 if (fPrealignment) delete fPrealignment; fPrealignment = 0;
170 if (fConstrRef) delete fConstrRef; fConstrRef = 0;
171 if (fCorrectSDD) delete fCorrectSDD; fCorrectSDD = 0;
172 if (fInitialRecSDD) delete fInitialRecSDD; fInitialRecSDD = 0;
173 fTrackBuff.Delete();
174 fConstraints.Delete();
7b85e477 175 fMilleModule.Delete();
176 fSuperModule.Delete();
177 if (--fgInstanceID==0) fgInstance = 0;
178}
179
180///////////////////////////////////////////////////////////////////////
18986853 181
182
6526a72c 183TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
184{
45993b23 185 // read new record from config file
6526a72c 186 TString record;
187 static TObjArray* recElems = 0;
188 if (recElems) {delete recElems; recElems = 0;}
189 //
190 TString keyws = recTitle;
191 if (!keyws.IsNull()) {
192 keyws.ToUpper();
193 // keyws += " ";
194 }
195 while (record.Gets(stream)) {
196 int cmt=record.Index("#");
197 if (cmt>=0) record.Remove(cmt); // skip comment
198 record.ReplaceAll("\t"," ");
199 record.ReplaceAll("\r"," ");
200 record.Remove(TString::kBoth,' ');
201 if (record.IsNull()) continue; // nothing to decode
202 if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
203 //
204 recElems = record.Tokenize(" ");
205 recTitle = recElems->At(0)->GetName();
206 recTitle.ToUpper();
207 recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
208 break;
209 }
210 if (rew || !recElems) rewind(stream);
211 return recElems;
212}
7b85e477 213
18986853 214//________________________________________________________________________________________________________
215Int_t AliITSAlignMille2::CheckConfigRecords(FILE* stream)
216{
217 TString record,recTitle;
218 int lineCnt = 0;
219 rewind(stream);
220 while (record.Gets(stream)) {
221 int cmt=record.Index("#");
222 lineCnt++;
223 if (cmt>=0) record.Remove(cmt); // skip comment
224 record.ReplaceAll("\t"," ");
225 record.ReplaceAll("\r"," ");
226 record.Remove(TString::kBoth,' ');
227 if (record.IsNull()) continue; // nothing to decode
228 // extract keyword
229 int spc = record.Index(" ");
230 if (spc>0) recTitle = record(0,spc);
231 else recTitle = record;
232 recTitle.ToUpper();
233 Bool_t strOK = kFALSE;
234 for (int ik=kNKeyWords;ik--;) if (recTitle == kRecKeys[ik]) {strOK = kTRUE; break;}
235 if (strOK) continue;
236 //
237 AliError(Form("Unknown keyword %s at line %d",
238 recTitle.Data(),lineCnt));
239 return -1;
240 //
241 }
242 //
243 rewind(stream);
244 return 0;
245}
246
247
6526a72c 248//________________________________________________________________________________________________________
249Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
45993b23 250{
251 // return 0 if success
252 // 1 if error in module index or voluid
6526a72c 253 //
7b85e477 254 FILE *pfc=fopen(cfile,"r");
255 if (!pfc) return -1;
6526a72c 256 //
257 TString record,recTitle,recOpt,recExt;
258 Int_t nrecElems,irec;
259 TObjArray *recArr=0;
260 //
261 fNModules = 0;
262 Bool_t stopped = kFALSE;
263 //
18986853 264 if (CheckConfigRecords(pfc)<0) return -1;
265 //
6526a72c 266 while(1) {
7b85e477 267 //
6526a72c 268 // ============= 1: we read some obligatory records in predefined order ================
269 //
18986853 270 recTitle = kRecKeys[kGeomFile];
6526a72c 271 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
272 (fGeometryFileName=recOpt).IsNull() ||
273 gSystem->AccessPathName(recOpt.Data()) ||
274 InitGeometry() )
275 { AliError("Failed to find/load Geometry"); stopped = kTRUE; break;}
7b85e477 276 //
18986853 277 recTitle = kRecKeys[kSuperModileFile];
6526a72c 278 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
279 recOpt.IsNull() ||
280 gSystem->AccessPathName(recOpt.Data()) ||
281 LoadSuperModuleFile(recOpt.Data()))
282 { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
7b85e477 283 //
18986853 284 recTitle = kRecKeys[kConstrRefFile]; // LOCAL_CONSTRAINTS are defined wrt these deltas
6526a72c 285 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
286 if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
287 else if (gSystem->AccessPathName(recOpt.Data()) || SetConstraintWrtRef(recOpt.Data()) )
288 { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
7b85e477 289 }
6526a72c 290 //
18986853 291 recTitle = kRecKeys[kPrealignFile];
6526a72c 292 if ( GetConfigRecord(pfc,recTitle,recOpt,1) )
293 if ( (fPreAlignmentFileName=recOpt).IsNull() ||
294 gSystem->AccessPathName(recOpt.Data()) ||
295 ApplyToGeometry())
296 { AliError(Form("Failed to load Prealignment file %s",recOpt.Data())); stopped = kTRUE; break;}
7b85e477 297 //
18986853 298 recTitle = kRecKeys[kPreCalSDDFile];
6526a72c 299 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
300 if ( recOpt.IsNull() || gSystem->AccessPathName(recOpt.Data()) ) {stopped = kTRUE; break;}
301 AliInfo(Form("Using %s for SDD precalibration",recOpt.Data()));
302 TFile* precfi = TFile::Open(recOpt.Data());
303 if (!precfi->IsOpen()) {stopped = kTRUE; break;}
304 fCorrectSDD = (AliITSresponseSDD*)precfi->Get("AliITSresponseSDD");
305 precfi->Close();
306 delete precfi;
307 if (!fCorrectSDD) {AliError("Precalibration SDD object is not found"); stopped = kTRUE; break;}
7b85e477 308 }
309 //
18986853 310 recTitle = kRecKeys[ kInitCalSDDFile ];
6526a72c 311 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
312 if ( recOpt.IsNull() || gSystem->AccessPathName(recOpt.Data()) ) {stopped = kTRUE; break;}
313 AliInfo(Form("Using %s as SDD calibration used in TrackPoints",recOpt.Data()));
314 TFile* precf = TFile::Open(recOpt.Data());
315 if (!precf->IsOpen()) {stopped = kTRUE; break;}
316 fInitialRecSDD = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
317 precf->Close();
318 delete precf;
319 if (!fInitialRecSDD) {AliError("Initial Calibration SDD object is not found"); stopped = kTRUE; break;}
7b85e477 320 }
321 //
18986853 322 recTitle = kRecKeys[ kGlobalDeltas ];
6526a72c 323 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
7b85e477 324 //
6526a72c 325 // =========== 2: see if there are local gaussian constraints defined =====================
326 // Note that they should be loaded before the modules declaration
7b85e477 327 //
18986853 328 recTitle = kRecKeys[ kConstrLocal ];
6526a72c 329 while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
330 nrecElems = recArr->GetLast()+1;
331 if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
332 if (GetConstraint(recOpt.Data())) {
333 AliError(Form("Existing constraint %s repeated",recOpt.Data()));
334 stopped = kTRUE; break;
335 }
336 recExt = recArr->At(2)->GetName();
337 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
338 double val = recExt.Atof();
339 recExt = recArr->At(3)->GetName();
340 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
341 double err = recExt.Atof();
342 int nwgh = nrecElems - 4;
343 double *wgh = new double[nwgh];
344 for (nwgh=0,irec=4;irec<nrecElems;irec++) {
345 recExt = recArr->At(irec)->GetName();
346 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
347 wgh[nwgh++] = recExt.Atof();
348 }
349 if (stopped) {delete[] wgh; break;}
350 //
351 ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
352 delete[] wgh;
353 //
354 } // end while for loop over local constraints
355 if (stopped) break;
7b85e477 356 //
6526a72c 357 // =========== 3: now read modules to align ===================================
7b85e477 358 //
6526a72c 359 rewind(pfc);
360 while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
18986853 361 if (!(recTitle==kRecKeys[ kModVolID ] || recTitle==kRecKeys[ kModIndex ])) continue;
6526a72c 362 // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ] extra params]
363 // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
364 // sig* is the scaling parameters for the errors of the clusters of this module
365 // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
366 //
367 nrecElems = recArr->GetLast()+1;
368 if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
369 int idx = recOpt.Atoi();
370 UShort_t voluid = (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
371 AliITSAlignMille2Module* mod = 0;
372 //
373 if (voluid>=kMinITSSupeModuleID) { // custom supermodule
374 for (int j=0; j<fNSuperModules; j++) {
375 if (voluid==GetSuperModule(j)->GetVolumeID()) {
376 mod = new AliITSAlignMille2Module(*GetSuperModule(j));
377 // the matrix might be updated in case some prealignment was applied, check
378 TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
379 if (mup) *(mod->GetMatrix()) = *mup;
380 fMilleModule.AddAtAndExpand(mod,fNModules);
381 break;
382 }
383 }
384 }
385 else if (idx<=kMaxITSSensVID) {
386 mod = new AliITSAlignMille2Module(voluid);
387 fMilleModule.AddAtAndExpand(mod,fNModules);
388 }
389 if (!mod) {stopped = kTRUE; break;} // bad volid
390 //
391 // geometry variation settings
392 for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
393 irec = i+2;
394 if (irec >= nrecElems) break;
395 recExt = recArr->At(irec)->GetName();
396 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
397 mod->SetFreeDOF(i, recExt.Atof() );
398 }
399 if (stopped) break;
400 //
401 // scaling factors for cluster errors
402 // first set default ones
403 for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);
404 for (int i=0;i<3;i++) {
405 irec = i+8;
406 if (irec >= nrecElems) break;
407 recExt = recArr->At(irec)->GetName();
408 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
409 mod->SetSigmaFactor(i, recExt.Atof() );
410 }
411 if (stopped) break;
412 //
413 mod->SetGeomParamsGlobal(fUseGlobalDelta);
414 // now comes special detectors treatment
415 if (mod->IsSDD()) {
416 double vl = 0;
417 if (nrecElems>11) {
418 recExt = recArr->At(11)->GetName();
419 if (recExt.IsFloat()) vl = recExt.Atof();
420 else {stopped = kTRUE; break;}
421 irec = 11;
422 }
423 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
424 //
425 vl = 0;
426 if (nrecElems>12) {
427 recExt = recArr->At(12)->GetName();
428 if (recExt.IsFloat()) vl = recExt.Atof();
429 else {stopped = kTRUE; break;}
430 irec = 12;
431 }
432 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFDV,vl);
433 }
434 //
435 mod->SetUniqueID(fNModules);
436 mod->EvaluateDOF();
437 fNModules++;
438 //
439 // now check if there are local constraints on this module
440 for (++irec;irec<nrecElems;irec++) {
441 recExt = recArr->At(irec)->GetName();
442 if (recExt.IsFloat()) {stopped=kTRUE;break;}
443 AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
444 if (!cstr) {
445 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
446 stopped=kTRUE;
447 break;
448 }
449 cstr->AddModule(mod);
450 }
451 if (stopped) break;
452 } // end while for loop over modules
453 if (stopped) break;
7b85e477 454 //
6526a72c 455 if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}
456 BuildHierarchy(); // preprocess loaded modules
7b85e477 457 //
6526a72c 458 // =========== 4: the rest may come in arbitrary order =======================================
459 rewind(pfc);
460 while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
461 //
462 nrecElems = recArr->GetLast()+1;
463 //
464 // some simple flags -----------------------------------------------------------------------
465 //
18986853 466 if (recTitle == kRecKeys[ kPseudoParents ]) SetAllowPseudoParents(kTRUE);
6526a72c 467 //
468 // some optional parameters ----------------------------------------------------------------
18986853 469 else if (recTitle == kRecKeys[ kTrackFitMethod ]) {
6526a72c 470 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
471 SetInitTrackParamsMeth(recOpt.Atoi());
7b85e477 472 }
6526a72c 473 //
18986853 474 else if (recTitle == kRecKeys[ kMinPntTrack ]) {
6526a72c 475 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
476 fMinNPtsPerTrack = recOpt.Atoi();
7b85e477 477 }
7b85e477 478 //
18986853 479 else if (recTitle == kRecKeys[ kNStDev ]) {
6526a72c 480 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
481 fNStdDev = (Int_t)recOpt.Atof();
482 }
7b85e477 483 //
18986853 484 else if (recTitle == kRecKeys[ kResCutInit ]) {
6526a72c 485 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
486 fResCutInitial = recOpt.Atof();
7b85e477 487 }
6526a72c 488 //
18986853 489 else if (recTitle == kRecKeys[ kResCutOther ]) {
6526a72c 490 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
491 fResCut = recOpt.Atof();
7b85e477 492 }
7b85e477 493 //
18986853 494 else if (recTitle == kRecKeys[ kLocalSigFactor ]) { //-------------------------
6526a72c 495 for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
496 fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
497 if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
498 }
499 if (stopped) break;
500 }
7b85e477 501 //
18986853 502 else if (recTitle == kRecKeys[ kStartFactor ]) { //-------------------------
6526a72c 503 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
504 fStartFac = recOpt.Atof();
505 }
506 //
18986853 507 else if (recTitle == kRecKeys[ kBField ]) { //-------------------------
6526a72c 508 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
509 fBField = recOpt.Atof();
510 if (fBField>0) {
511 fBOn = kTRUE;
512 fNLocal = 5; // helices
513 fRieman = new AliTrackFitterRieman();
514 }
515 else {
516 fBField = 0.0;
517 fBOn = kFALSE;
518 fNLocal = 4;
519 }
520 }
521 //
18986853 522 else if (recTitle == kRecKeys[ kSparseMatrix ]) { // matrix solver type
6526a72c 523 //
524 AliMillePede2::SetGlobalMatSparse(kTRUE);
525 if (recOpt.IsNull()) continue;
526 // solver type and settings
527 if (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
528 else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
529 else {stopped = kTRUE; break;}
530 //
531 if (nrecElems>=3) { // preconditioner type
532 recExt = recArr->At(2)->GetName();
533 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
534 AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
535 }
536 //
537 if (nrecElems>=4) { // tolerance
538 recExt = recArr->At(3)->GetName();
539 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
540 AliMillePede2::SetMinResTol( recExt.Atof() );
541 }
542 //
543 if (nrecElems>=5) { // maxIter
544 recExt = recArr->At(4)->GetName();
545 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
546 AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
547 }
548 }
549 //
18986853 550 else if (recTitle == kRecKeys[ kRequirePoint ]) { //-------------------------
6526a72c 551 // syntax: REQUIRE_POINT where ndet updw nreqpts
552 // where = LAYER or DETECTOR
553 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
554 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
555 // nreqpts = minimum number of points of that type
556 if (nrecElems>=5) {
557 recOpt.ToUpper();
558 int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
559 int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
560 int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
561 fRequirePoints = kTRUE;
562 if (recOpt == "LAYER") {
563 if (lr<0 || lr>5) {stopped = kTRUE; break;}
564 if (hb>0) fNReqLayUp[lr] = np;
565 else if (hb<0) fNReqLayDown[lr] = np;
566 else fNReqLay[lr] = np;
567 }
568 else if (recOpt == "DETECTOR") {
569 if (lr<0 || lr>2) {stopped = kTRUE; break;}
570 if (hb>0) fNReqDetUp[lr] = np;
571 else if (hb<0) fNReqDetDown[lr] = np;
572 else fNReqDet[lr] = np;
573 }
574 else {stopped = kTRUE; break;}
575 }
576 else {stopped = kTRUE; break;}
577 }
578 //
579 // global constraints on the subunits/orphans
18986853 580 else if (recTitle == kRecKeys[ kConstrOrphans ]) { //------------------------
6526a72c 581 // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
582 if (nrecElems<4) {stopped = kTRUE; break;}
583 recExt = recArr->At(2)->GetName();
584 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
585 double val = recExt.Atof();
586 UInt_t pattern = 0;
587 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
588 recExt = recArr->At(irec)->GetName();
589 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
590 pattern |= 0x1 << recExt.Atoi();
591 }
592 if (stopped) break;
593 if (recOpt == "MEAN") ConstrainOrphansMean(val,pattern);
594 else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
595 else {stopped = kTRUE; break;}
596 }
597 //
18986853 598 else if (recTitle == kRecKeys[ kConstrSubunits ]) { //------------------------
6526a72c 599 // expect ONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
600 if (nrecElems<5) {stopped = kTRUE; break;}
601 recExt = recArr->At(2)->GetName();
602 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
603 double val = recExt.Atof();
604 UInt_t pattern = 0;
605 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
606 recExt = recArr->At(irec)->GetName();
607 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
608 int parid = recExt.Atoi();
609 if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
610 else break; // list of params is over
611 }
612 if (stopped) break;
613 //
614 Bool_t meanC;
615 if (recOpt == "MEAN") meanC = kTRUE;
616 else if (recOpt == "MEDIAN") meanC = kFALSE;
617 else {stopped = kTRUE; break;}
618 //
619 int curID = -1;
620 int rangeStart = -1;
621 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
622 recExt = recArr->At(irec)->GetName();
623 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
624 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
625 else curID = recExt.Atoi();
626 //
627 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
628 // this was a range start or single
629 int start;
630 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
631 else start = curID; // create constraint either for single module (or 1st in the range)
632 for (int id=start;id<=curID;id++) {
633 int id0 = IsVIDDefined(id);
634 if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
635 if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
636 else ConstrainModuleSubUnitsMedian(id0,val,pattern);
637 }
638 }
639 if (rangeStart>=0) stopped = kTRUE; // unfinished range
640 if (stopped) break;
641 }
642 //
643 // association of modules with local constraints
18986853 644 else if (recTitle == kRecKeys[ kApplyConstr ]) { //------------------------
6526a72c 645 // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
646 if (nrecElems<3) {stopped = kTRUE; break;}
647 int nmID0=-1,nmID1=-1;
648 for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
649 recExt = recArr->At(irec)->GetName();
650 if (recExt.IsFloat()) break;
651 // check if such a constraint was declared
652 if (!GetConstraint(recExt.Data())) {
653 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
654 stopped=kTRUE;
655 break;
656 }
657 if (nmID0<0) nmID0 = irec;
658 nmID1 = irec;
659 }
660 if (stopped) break;
661 //
662 if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
663 //
664 // now read the list of modules to constrain
665 int curID = -1;
666 int rangeStart = -1;
667 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
668 recExt = recArr->At(irec)->GetName();
669 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
670 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
671 else curID = recExt.Atoi();
672 //
673 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
674 //
675 // this was a range start or single
676 int start;
677 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
678 else start = curID; // create constraint either for single module (or 1st in the range)
679 for (int id=start;id<=curID;id++) {
680 AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
681 if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
682 for (int nmid=nmID0;nmid<=nmID1;nmid++)
683 ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
684 }
685 }
686 if (rangeStart>=0) stopped = kTRUE; // unfinished range
687 if (stopped) break;
688 }
689 //
690 else continue; // already processed record
691 //
692 } // end of while loop 4 over the various params
7b85e477 693 //
6526a72c 694 break;
695 } // end of while(1) loop
7b85e477 696 //
697 fclose(pfc);
6526a72c 698 if (stopped) {
699 AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
700 return -1;
701 }
702 //
703 fIsConfigured = kTRUE;
704 return 0;
705}
706
707//________________________________________________________________________________________________________
708void AliITSAlignMille2::BuildHierarchy()
709{
710 // build the hieararhy of the modules to align
7b85e477 711 //
6526a72c 712 if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
713 AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
714 "Since Deltas are local, switching to NoPseudoParents");
715 SetAllowPseudoParents(kFALSE);
716 }
7b85e477 717 // set parent/child relationship for modules to align
6526a72c 718 AliInfo("Setting parent/child relationships\n");
7b85e477 719 //
6526a72c 720 // 1) child -> parent reference
721 for (int ipar=0;ipar<fNModules;ipar++) {
7b85e477 722 AliITSAlignMille2Module* parent = GetMilleModule(ipar);
6526a72c 723 if (parent->IsSensor()) continue; // sensor cannot be a parent
7b85e477 724 //
6526a72c 725 for (int icld=0;icld<fNModules;icld++) {
7b85e477 726 if (icld==ipar) continue;
727 AliITSAlignMille2Module* child = GetMilleModule(icld);
728 if (!child->BelongsTo(parent)) continue;
6526a72c 729 // child cannot have more sensors than the parent
730 if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
7b85e477 731 //
732 AliITSAlignMille2Module* parOld = child->GetParent();
6526a72c 733 // is this parent candidate closer than the old parent ?
7b85e477 734 if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
735 child->SetParent(parent);
736 }
737 //
738 }
739 //
6526a72c 740 // add parent -> children reference
741 for (int icld=0;icld<fNModules;icld++) {
742 AliITSAlignMille2Module* child = GetMilleModule(icld);
743 AliITSAlignMille2Module* parent = child->GetParent();
744 if (parent) parent->AddChild(child);
745 }
746 //
7b85e477 747 // reorder the modules in such a way that parents come first
6526a72c 748 for (int icld=0;icld<fNModules;icld++) {
057304f3 749 AliITSAlignMille2Module* child = GetMilleModule(icld);
750 AliITSAlignMille2Module* parent;
6526a72c 751 while ( (parent=child->GetParent()) && (parent->GetUniqueID()>child->GetUniqueID()) ) {
7b85e477 752 // swap
753 fMilleModule[icld] = parent;
754 fMilleModule[parent->GetUniqueID()] = child;
755 child->SetUniqueID(parent->GetUniqueID());
756 parent->SetUniqueID(icld);
057304f3 757 child = parent;
7b85e477 758 }
057304f3 759 //
6526a72c 760 }
7b85e477 761 //
6526a72c 762 // Go over the child->parent chain and mark modules with explicitly provided sensors.
763 // If the sensors of the unit are explicitly declared, all undeclared sensors are
764 // suppresed in this unit.
765 for (int icld=fNModules;icld--;) {
7b85e477 766 AliITSAlignMille2Module* child = GetMilleModule(icld);
767 AliITSAlignMille2Module* parent = child->GetParent();
768 if (!parent) continue;
6526a72c 769 //
770 // check if this parent was already processed
771 if (!parent->AreSensorsProvided()) {
772 parent->DelSensitiveVolumes();
773 parent->SetSensorsProvided(kTRUE);
774 }
775 // reattach sensors to parent
776 for (int isc=child->GetNSensitiveVolumes();isc--;) {
777 UShort_t senVID = child->GetSensVolVolumeID(isc);
778 if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
7b85e477 779 }
780 }
781 //
7b85e477 782}
783
6526a72c 784// pepo
785//________________________________________________________________________________________________________
786void AliITSAlignMille2::SetCurrentModule(Int_t id)
787{
45993b23 788 // set the current supermodule
6526a72c 789 // new meaning
790 if (fMilleVersion>=2) {
791 fCurrentModule = GetMilleModule(id);
792 return;
793 }
794 // old meaning
795 if (fMilleVersion<=1) {
796 Int_t index=id;
797 /// set as current the SuperModule that contains the 'index' sens.vol.
798 if (index<0 || index>2197) {
799 AliInfo("index does not correspond to a sensitive volume!");
800 return;
801 }
802 UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
803 Int_t k=IsContained(voluid);
804 if (k>=0){
805 fCluster.SetVolumeID(voluid);
806 fCluster.SetXYZ(0,0,0);
807 InitModuleParams();
808 }
809 else
810 AliInfo(Form("module %d not defined\n",index));
811 }
812}
813// endpepo
814//________________________________________________________________________________________________________
7b85e477 815void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts)
816{
817 // set minimum number of points in specific detector or layer
818 // where = LAYER or DETECTOR
819 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
820 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
821 // nreqpts = minimum number of points of that type
822 ndet--;
823 if (strstr(where,"LAYER")) {
824 if (ndet<0 || ndet>5) return;
825 if (updw>0) fNReqLayUp[ndet]=nreqpts;
826 else if (updw<0) fNReqLayDown[ndet]=nreqpts;
827 else fNReqLay[ndet]=nreqpts;
828 fRequirePoints=kTRUE;
829 }
830 else if (strstr(where,"DETECTOR")) {
831 if (ndet<0 || ndet>2) return;
832 if (updw>0) fNReqDetUp[ndet]=nreqpts;
833 else if (updw<0) fNReqDetDown[ndet]=nreqpts;
834 else fNReqDet[ndet]=nreqpts;
835 fRequirePoints=kTRUE;
836 }
837}
838
6526a72c 839//________________________________________________________________________________________________________
840Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname)
841{
7b85e477 842 /// index from symname
843 if (!symname) return -1;
844 for (Int_t i=0;i<=kMaxITSSensID; i++) {
845 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
846 }
847 return -1;
848}
849
6526a72c 850//________________________________________________________________________________________________________
851Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid)
852{
7b85e477 853 /// index from volume ID
854 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
855 if (lay<1|| lay>6) return -1;
856 Int_t idx=Int_t(voluid)-2048*lay;
857 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
858 for (Int_t ilay=1; ilay<lay; ilay++)
859 idx += AliGeomManager::LayerSize(ilay);
860 return idx;
861}
862
6526a72c 863//________________________________________________________________________________________________________
864UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname)
865{
7b85e477 866 /// volume ID from symname
867 /// works for sensitive volumes only
868 if (!symname) return 0;
869
870 for (UShort_t voluid=2000; voluid<13300; voluid++) {
871 Int_t modId;
872 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
873 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
874 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
875 }
876 }
877
878 return 0;
879}
880
6526a72c 881//________________________________________________________________________________________________________
882UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index)
883{
7b85e477 884 /// volume ID from index
885 if (index<0) return 0;
886 if (index<2198)
887 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
888 else {
889 for (int i=0; i<fNSuperModules; i++) {
890 if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
891 }
892 }
893 return 0;
894}
895
6526a72c 896//________________________________________________________________________________________________________
897Int_t AliITSAlignMille2::InitGeometry()
898{
7b85e477 899 /// initialize geometry
6526a72c 900 AliInfo("Loading initial geometry");
7b85e477 901 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
902 fGeoManager = AliGeomManager::GetGeometry();
903 if (!fGeoManager) {
904 AliInfo("Couldn't initialize geometry");
6526a72c 905 return -1;
906 }
907 return 0;
908}
909
910//________________________________________________________________________________________________________
911Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname)
912{
913 // Load the global deltas from this file. The local gaussian constraints on some modules
914 // will be defined with respect to the deltas from this reference file, converted to local
915 // delta format. Note: conversion to local format requires reloading the geometry!
916 //
917 AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
918 if (!fGeoManager) return -1;
919 fConstrRefFileName = reffname;
920 if (fConstrRefFileName == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
921 fConstrRef = new TClonesArray("AliAlignObjParams",1);
922 return 0;
923 }
924 TFile *pref = TFile::Open(fConstrRefFileName.Data());
925 if (!pref->IsOpen()) return -2;
926 fConstrRef = (TClonesArray*)pref->Get("ITSAlignObjs");
927 pref->Close();
928 delete pref;
929 if (!fConstrRef) {
930 AliError(Form("Did not find reference prealignment deltas in %s",reffname));
931 return -1;
932 }
933 //
934 // we need ideal geometry to convert global deltas to local ones
935 if (fUsePreAlignment) {
936 AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
937 return -1;
938 }
939 //
940 AliInfo("Converting global reference deltas to local ones");
941 Int_t nprea = fConstrRef->GetEntriesFast();
942 for (int ix=0; ix<nprea; ix++) {
943 AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
944 if (!preo->ApplyToGeometry()) return -1;
7b85e477 945 }
6526a72c 946 //
947 // now convert the global reference deltas to local ones
948 for (int i=fConstrRef->GetEntriesFast();i--;) {
949 AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
950 TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
951 if (!mupd) { // this is not alignable entry, need to look in the supermodules
952 for (int im=fNSuperModules;im--;) {
953 AliITSAlignMille2Module* mod = GetSuperModule(im);
954 if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
955 mupd = mod->GetMatrix();
956 break;
957 }
958 if (!mupd) {
959 AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
960 return -1;
961 }
962 }
963 TGeoHMatrix preMat;
964 preo->GetMatrix(preMat); // Delta_Glob
965 TGeoHMatrix tmpMat = *mupd; // Delta_Glob * Delta_Glob_Par * M
966 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
967 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
968 preo->SetMatrix(tmpMat); // local corrections
969 }
970 //
971 // we need to reload the geometry spoiled by this reference deltas...
972 delete fGeoManager;
973 AliInfo("Reloading initial geometry");
974 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
975 fGeoManager = AliGeomManager::GetGeometry();
976 return 0;
7b85e477 977}
978
6526a72c 979//________________________________________________________________________________________________________
980void AliITSAlignMille2::Init()
7b85e477 981{
45993b23 982 // perform global initialization
6526a72c 983 //
984 if (fIsMilleInit) {
985 AliInfo("Millepede has been already initialized!");
986 return;
987 }
988 // range constraints in such a way that the childs are constrained before their parents
989 // orphan constraints come last
990 for (int ic=0;ic<GetNConstraints();ic++) {
991 for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
992 AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
993 AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
994 if (cst0->GetModuleID()<cst1->GetModuleID()) {
995 // swap
996 fConstraints[ic] = cst1;
997 fConstraints[ic1] = cst0;
998 }
999 }
1000 }
1001 //
1002 if (!GetUseGlobalDelta()) {
1003 AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
1004 for (int imd=fNModules;imd--;) {
1005 AliITSAlignMille2Module* mod = GetMilleModule(imd);
1006 int npar = mod->GetNParTot();
1007 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1008 for (int ipar=0;ipar<npar;ipar++) {
1009 if (!mod->IsFreeDOF(ipar)) continue;
1010 mod->SetParOffset(ipar,fNGlobal++);
1011 }
1012 }
1013 }
1014 else {
1015 // init millepede, decide which parameters are to be fitted explicitly
1016 for (int imd=fNModules;imd--;) {
1017 AliITSAlignMille2Module* mod = GetMilleModule(imd);
1018 int npar = mod->GetNParTot();
1019 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1020 for (int ipar=0;ipar<npar;ipar++) {
1021 if (!mod->IsFreeDOF(ipar)) continue; // fixed
1022 //
1023 int nFreeInstances = 0;
1024 //
1025 AliITSAlignMille2Module* parent = mod;
1026 Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
1027 //
45993b23 1028 Bool_t addToFit = kFALSE;
6526a72c 1029 // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
1030 // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
1031 // 2) the same applies to all of its parents
1032 // 3) it has at least 1 unconstrained direct child
1033 while(parent) {
1034 if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
1035 nFreeInstances++;
1036 if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
45993b23 1037 if (cstGauss) addToFit = kTRUE;
6526a72c 1038 parent = parent->GetParent();
1039 }
1040 if (nFreeInstances>1) {
1041 AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
1042 "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
1043 exit(1);
1044 }
1045 //
1046 // i) Are PseudoParents allowed?
45993b23 1047 if (!PseudoParentsAllowed()) addToFit = kTRUE;
6526a72c 1048 // ii) check if this module has no child with such a free parameter. Since the order of this check
1049 // goes from child to parent, by this moment such a parameter must have been already added
45993b23 1050 else if (!IsParModFamilyVaried(mod,ipar)) addToFit = kTRUE; // no varied children at all
1051 else if (!IsParFamilyFree(mod,ipar,1)) addToFit = kTRUE; // no unconstrained direct children
6526a72c 1052 // otherwise the value of this parameter can be extracted from simple contraint and the values of
1053 // the relevant parameters of its children the fit is done. Hence it is not included
45993b23 1054 if (!addToFit) continue;
6526a72c 1055 //
1056 // shall add this parameter to explicit fit
1057 // printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
1058 mod->SetParOffset(ipar,fNGlobal++);
1059 }
1060 }
1061 }
1062 //
1063 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
1064 fGlobalDerivatives = new Double_t[fNGlobal];
1065 memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1066 //
1067 fMillepede->InitMille(fNGlobal,fNLocal,fNStdDev,fResCut,fResCutInitial);
7b85e477 1068 fIsMilleInit = kTRUE;
6526a72c 1069 //
1070 ResetLocalEquation();
1071 AliInfo("Parameters initialized to zero");
1072 //
7b85e477 1073 /// Fix non free parameters
1074 for (Int_t i=0; i<fNModules; i++) {
6526a72c 1075 AliITSAlignMille2Module* mod = GetMilleModule(i);
1076 for (Int_t j=0; j<mod->GetNParTot(); j++) {
1077 if (mod->GetParOffset(j)<0) continue; // not varied
1078 FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1079 fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1080 }
7b85e477 1081 }
1082 //
1083 // Set iterations
6526a72c 1084 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
1085 //
7b85e477 1086}
1087
6526a72c 1088//________________________________________________________________________________________________________
1089void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma)
1090{
7b85e477 1091 /// Constrain equation defined by par to value
6526a72c 1092 if (!fIsMilleInit) Init();
1093 fMillepede->SetGlobalConstraint(par, value, sigma);
7b85e477 1094 AliInfo("Adding constraint");
1095}
1096
6526a72c 1097//________________________________________________________________________________________________________
1098void AliITSAlignMille2::InitGlobalParameters(Double_t *par)
1099{
7b85e477 1100 /// Initialize global parameters with par array
6526a72c 1101 if (!fIsMilleInit) Init();
7b85e477 1102 fMillepede->SetGlobalParameters(par);
1103 AliInfo("Init Global Parameters");
1104}
6526a72c 1105
1106//________________________________________________________________________________________________________
1107void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value)
1108{
7b85e477 1109 /// Parameter iPar is encourage to vary in [-value;value].
1110 /// If value == 0, parameter is fixed
1111 if (!fIsMilleInit) {
1112 AliInfo("Millepede has not been initialized!");
1113 return;
1114 }
1115 fMillepede->SetParSigma(iPar, value);
1116 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
1117}
1118
6526a72c 1119//________________________________________________________________________________________________________
7b85e477 1120void AliITSAlignMille2::ResetLocalEquation()
1121{
1122 /// Reset the derivative vectors
1123 for(int i=fNLocal;i--;) fLocalDerivatives[i] = 0.0;
6526a72c 1124 memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
7b85e477 1125}
1126
6526a72c 1127//________________________________________________________________________________________________________
7b85e477 1128Int_t AliITSAlignMille2::ApplyToGeometry()
1129{
1130 // apply starting realignment to ideal geometry
6526a72c 1131 AliInfo(Form("Using %s for prealignment",fPreAlignmentFileName.Data()));
7b85e477 1132 if (!fGeoManager) return -1;
6526a72c 1133 TFile *pref = TFile::Open(fPreAlignmentFileName.Data());
7b85e477 1134 if (!pref->IsOpen()) return -2;
1135 fPrealignment = (TClonesArray*)pref->Get("ITSAlignObjs");
1136 if (!fPrealignment) return -3;
1137 Int_t nprea = fPrealignment->GetEntriesFast();
1138 AliInfo(Form("Array of input misalignments with %d entries",nprea));
1139 //
1140 for (int ix=0; ix<nprea; ix++) {
1141 AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1142 Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1143 if (index>=0) {
1144 if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1145 fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1146 }
1147 //TString nms = preo->GetSymName();
1148 //if (!nms.Contains("Ladder")) continue; //RRR
6526a72c 1149 //printf("Applying#%4d %s\n",ix,preo->GetSymName());
7b85e477 1150 if (!preo->ApplyToGeometry()) return -4;
1151 }
1152 //
1153 pref->Close();
1154 delete pref;
1155 //
1156 fUsePreAlignment = kTRUE;
1157 return 0;
1158}
1159
6526a72c 1160//________________________________________________________________________________________________________
7b85e477 1161Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1162{
45993b23 1163 // quality factors from prealignment
7b85e477 1164 if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1165 return fPreAlignQF[index]-1;
1166}
1167
6526a72c 1168//________________________________________________________________________________________________________
1169AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp)
1170{
7b85e477 1171 /// create a new AliTrackPointArray keeping only defined modules
1172 /// move points according to a given prealignment, if any
1173 /// sort alitrackpoints w.r.t. global Y direction, if selected
6526a72c 1174 const double kTiny = 1E-12;
1175 //
7b85e477 1176 AliTrackPointArray *atps=NULL;
1177 Int_t idx[20];
1178 Int_t npts=atp->GetNPoints();
1179
1180 /// checks if AliTrackPoints belong to defined modules
1181 Int_t ngoodpts=0;
1182 Int_t intidx[20];
1183
1184 for (int j=0; j<npts; j++) {
6526a72c 1185 intidx[j] = IsVIDContained(atp->GetVolumeID()[j]);
7b85e477 1186 if (intidx[j]>=0) ngoodpts++;
1187 }
6526a72c 1188 AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
7b85e477 1189
1190 // reject track if not enough points are left
1191 if (ngoodpts<fMinNPtsPerTrack) {
1192 AliInfo("Track with not enough points!");
1193 return NULL;
1194 }
1195 // >> RS
1196 AliTrackPoint p;
1197 // check points in specific places
1198 if (fRequirePoints) {
1199 Int_t nlayup[6],nlaydown[6],nlay[6];
1200 Int_t ndetup[3],ndetdown[3],ndet[3];
1201 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1202 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1203
1204 for (int i=0; i<npts; i++) {
1205 // skip not defined points
1206 if (intidx[i]<0) continue;
1207 Float_t xx=atp->GetX()[i];
1208 Float_t yy=atp->GetY()[i];
1209 Float_t r=TMath::Sqrt(xx*xx + yy*yy);
1210 int lay=-1;
1211 if (r<5) lay=0;
1212 else if (r>5 && r<10) lay=1;
1213 else if (r>10 && r<18) lay=2;
1214 else if (r>18 && r<30) lay=3;
1215 else if (r>30 && r<40) lay=4;
1216 else if (r>40) lay=5;
1217 if (lay<0) continue;
1218 int det=lay/2;
1219 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
1220
1221 if (yy>=0.0) { // UP point
1222 nlayup[lay]++;
1223 nlay[lay]++;
1224 ndetup[det]++;
1225 ndet[det]++;
1226 }
1227 else {
1228 nlaydown[lay]++;
1229 nlay[lay]++;
1230 ndetdown[det]++;
1231 ndet[det]++;
1232 }
1233 }
1234
1235 // checks minimum values
1236 Bool_t isok=kTRUE;
1237 for (Int_t j=0; j<6; j++) {
1238 if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE;
1239 if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE;
1240 if (nlay[j]<fNReqLay[j]) isok=kFALSE;
1241 }
1242 for (Int_t j=0; j<3; j++) {
1243 if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE;
1244 if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE;
1245 if (ndet[j]<fNReqDet[j]) isok=kFALSE;
1246 }
1247 if (!isok) {
1248 AliDebug(2,Form("Track does not meet all location point requirements!"));
1249 return NULL;
1250 }
1251 }
7b85e477 1252 // build a new track with (sorted) (prealigned) good points
6526a72c 1253 atps = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
1254 if (!atps) {
1255 atps = new AliTrackPointArray(ngoodpts);
1256 fTrackBuff.AddAtAndExpand(atps,ngoodpts-fMinNPtsPerTrack);
1257 }
1258 //
1259 //
7b85e477 1260 for (int i=0; i<npts; i++) idx[i]=i;
1261 // sort track if required
6526a72c 1262 TMath::Sort(npts,atp->GetY(),idx); // sort descending...
1263 //
7b85e477 1264 Int_t npto=0;
1265 for (int i=0; i<npts; i++) {
1266 // skip not defined points
1267 if (intidx[idx[i]]<0) continue;
1268 atp->GetPoint(p,idx[i]);
1269
1270 // prealign point if required
1271 // get IDEAL matrix
6526a72c 1272 AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1273 TGeoHMatrix *svOrigMatrix = mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
7b85e477 1274 // get back real local coordinates: use OriginalGlobalMatrix because AliTrackPoints were written
1275 // with idel geometry
1276 Double_t pg[3],pl[3];
1277 pg[0]=p.GetX();
1278 pg[1]=p.GetY();
1279 pg[2]=p.GetZ();
1280 // printf("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1281 AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1282 svOrigMatrix->MasterToLocal(pg,pl);
1283
1284 AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]));
6526a72c 1285 //
1286 // this is a temporary code to extract the drift speed used for given point
1287 if (p.GetDriftTime()>0) { // RRR
1288 // calculate the drift speed
1289 int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());// - kSDDoffsID;
1290 fDriftTime0[npto] = fInitialRecSDD ? fInitialRecSDD->GetTimeZero(sid) : 0.;
1291 /*
1292 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(p.GetVolumeID());
1293 if (lay==3) fDriftTime0[npto] = pg[2]<0 ? 169.5 : 140.1;
1294 else if (lay==4) fDriftTime0[npto] = pg[2]<0 ? 158.3 : 139.0;
1295 else {
1296 AliError(Form("Strange layer %d for moduleID %d",lay,p.GetVolumeID()));
1297 exit(1);
1298 }
1299 */
1300 double tdif = p.GetDriftTime() - fDriftTime0[npto];
1301 if (tdif<=0) tdif = 1;
1302 double vdrift = (3.5085-TMath::Abs(pl[0]))/tdif;
1303 if (vdrift<0) vdrift = 0;
1304 //
1305 // TEMPORARY CORRECTION (if provided) -------------->>>
1306 if (fCorrectSDD) {
1307 float t0Upd = fCorrectSDD->GetTimeZero(sid);
1308 vdrift += fCorrectSDD->GetDeltaVDrift(sid);
1309 tdif = p.GetDriftTime() - t0Upd;
1310 // correct Xlocal
1311 pl[0] = TMath::Sign(3.5085 - vdrift*tdif,pl[0]);
1312 fDriftTime0[npto] = t0Upd;
1313 }
1314 // TEMPORARY CORRECTION (if provided) --------------<<<
1315 fDriftSpeed[npto] = TMath::Sign(vdrift,pl[0]);
1316 //
1317 /*
1318 printf("%d %+6.2f %+6.2f %+6.2f %+5.2f %+5.2f %+5.2f %+6.1f %+6.1f %+f %+f\n",
1319 p.GetVolumeID(),pg[0],pg[1],pg[2],pl[0],pl[1],pl[2],p.GetDriftTime(), fDriftTime0[npto], fDriftSpeed[npto],tdif);
1320 */
1321 }
7b85e477 1322
1323 // update covariance matrix
1324 TGeoHMatrix hcov;
1325 Double_t hcovel[9];
1326 hcovel[0]=double(p.GetCov()[0]);
1327 hcovel[1]=double(p.GetCov()[1]);
1328 hcovel[2]=double(p.GetCov()[2]);
1329 hcovel[3]=double(p.GetCov()[1]);
1330 hcovel[4]=double(p.GetCov()[3]);
1331 hcovel[5]=double(p.GetCov()[4]);
1332 hcovel[6]=double(p.GetCov()[2]);
1333 hcovel[7]=double(p.GetCov()[4]);
1334 hcovel[8]=double(p.GetCov()[5]);
1335 hcov.SetRotation(hcovel);
1336 // now rotate in local system
6526a72c 1337 // printf("\nErrMatGlob: before\n"); hcov.Print(""); //RRR
7b85e477 1338 hcov.Multiply(svOrigMatrix);
1339 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1340 // now hcov is LOCAL COVARIANCE MATRIX
6526a72c 1341 // apply sigma scaling
1342 // printf("\nErrMatLoc: before\n"); hcov.Print(""); //RRR
1343 Double_t *hcovscl = hcov.GetRotationMatrix();
1344 // for (int ir=3;ir--;) for (int ic=3;ic--;) hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1345 // RS TEMPORARY: nullify non-diagonal elements and sigY
1346 hcovscl[5] = 0;
1347 for (int ir=3;ir--;) for (int ic=3;ic--;) {
1348 if (ir==ic) {
1349 if (TMath::Abs(hcovscl[ir*3+ic])<kTiny) hcovscl[ir*3+ic] = 0.;
1350 else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1351 }
1352 else hcovscl[ir*3+ic] = 0;
1353 }
1354 //
1355 // printf("\nErrMatLoc: after\n"); hcov.Print(""); //RRR
1356 //
7b85e477 1357 if (fBug==1) {
1358 // correzione bug LAYER 5 SSD temporanea..
1359 int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1360 if (ssdidx>=500 && ssdidx<1248) {
1361 int ladder=(ssdidx-500)%22;
1362 if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1363 if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1364 }
1365 }
7b85e477 1366 /// get (evenctually prealigned) matrix of sens. vol.
6526a72c 1367 TGeoHMatrix *svMatrix = mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
7b85e477 1368 // modify global coordinates according with pre-aligment
1369 svMatrix->LocalToMaster(pl,pg);
1370 // now rotate in local system
1371 hcov.Multiply(&svMatrix->Inverse());
1372 hcov.MultiplyLeft(svMatrix);
1373 // hcov is back in GLOBAL RF
6526a72c 1374 // cure once more
1375 for (int ir=3;ir--;) for (int ic=3;ic--;) if (TMath::Abs(hcovscl[ir*3+ic])<kTiny) hcovscl[ir*3+ic] = 0.;
1376 // printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1377 //
7b85e477 1378 Float_t pcov[6];
6526a72c 1379 pcov[0]=hcovscl[0];
1380 pcov[1]=hcovscl[1];
1381 pcov[2]=hcovscl[2];
1382 pcov[3]=hcovscl[4];
1383 pcov[4]=hcovscl[5];
1384 pcov[5]=hcovscl[8];
7b85e477 1385
1386 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1387 // printf("New Gl coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1388 AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1389 atps->AddPoint(npto,&p);
6526a72c 1390 AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,atps->GetX()[npto],
1391 atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] ));
1392 // printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY());
7b85e477 1393 npto++;
1394 }
1395
1396 return atps;
1397}
1398
6526a72c 1399//________________________________________________________________________________________________________
1400AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp)
1401{
7b85e477 1402 /// sort alitrackpoints w.r.t. global Y direction
1403 AliTrackPointArray *atps=NULL;
1404 Int_t idx[20];
1405 Int_t npts=atp->GetNPoints();
1406 AliTrackPoint p;
1407 atps=new AliTrackPointArray(npts);
1408
1409 TMath::Sort(npts,atp->GetY(),idx);
1410
1411 for (int i=0; i<npts; i++) {
1412 atp->GetPoint(p,idx[i]);
1413 atps->AddPoint(i,&p);
1414 AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
1415 }
1416 return atps;
1417}
1418
6526a72c 1419//________________________________________________________________________________________________________
1420Int_t AliITSAlignMille2::GetCurrentLayer() const
1421{
45993b23 1422 // get current layer id
6526a72c 1423 if (!fGeoManager) {
1424 AliInfo("ITS geometry not initialized!");
1425 return -1;
1426 }
1427 return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1428}
7b85e477 1429
6526a72c 1430//________________________________________________________________________________________________________
1431Int_t AliITSAlignMille2::InitModuleParams()
1432{
7b85e477 1433 /// initialize geometry parameters for a given detector
1434 /// for current cluster (fCluster)
1435 /// fGlobalInitParam[] is set as:
1436 /// [tx,ty,tz,psi,theta,phi]
1437 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1438 /// *** At the moment: using Raffalele's angles definition ***
1439 ///
1440 /// return 0 if success
6526a72c 1441 /// If module is found but has no parameters to vary, return 1
7b85e477 1442
1443 if (!fGeoManager) {
1444 AliInfo("ITS geometry not initialized!");
1445 return -1;
1446 }
1447
1448 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1449
1450 // set the internal index (index in module list)
1451 UShort_t voluid=fCluster.GetVolumeID();
1452 //
1453 // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1454 Int_t k=fNModules-1;
1455 fCurrentModule = 0;
6526a72c 1456 // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules
1457 while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
7b85e477 1458 if (k<0) return -3;
6526a72c 1459 //
1460 /*
1461 // Check if the module has free params. If not, go over the parents
1462 AliITSAlignMille2Module* mdtmp = fCurrentModule;
1463 while (mdtmp && mdtmp->GetNParFree()==0) mdtmp = mdtmp->GetParent();
1464 if (!mdtmp) return 1; // nothing to vary here
1465 fCurrentModule = mdtmp;
1466 */
7b85e477 1467 //
1468 fModuleInitParam[0] = 0.0;
1469 fModuleInitParam[1] = 0.0;
1470 fModuleInitParam[2] = 0.0;
1471 fModuleInitParam[3] = 0.0; // psi (X)
1472 fModuleInitParam[4] = 0.0; // theta (Y)
1473 fModuleInitParam[5] = 0.0; // phi (Z)
6526a72c 1474 fModuleInitParam[6] = 0.0;
1475 fModuleInitParam[7] = 0.0;
7b85e477 1476 /// get (evenctually prealigned) matrix of sens. vol.
1477 TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1478
1479 fMeasGlo[0] = fCluster.GetX();
1480 fMeasGlo[1] = fCluster.GetY();
1481 fMeasGlo[2] = fCluster.GetZ();
1482 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1483 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
1484
1485 // set stdev from cluster
1486 TGeoHMatrix hcov;
1487 Double_t hcovel[9];
1488 hcovel[0]=double(fCluster.GetCov()[0]);
1489 hcovel[1]=double(fCluster.GetCov()[1]);
1490 hcovel[2]=double(fCluster.GetCov()[2]);
1491 hcovel[3]=double(fCluster.GetCov()[1]);
1492 hcovel[4]=double(fCluster.GetCov()[3]);
1493 hcovel[5]=double(fCluster.GetCov()[4]);
1494 hcovel[6]=double(fCluster.GetCov()[2]);
1495 hcovel[7]=double(fCluster.GetCov()[4]);
1496 hcovel[8]=double(fCluster.GetCov()[5]);
1497 hcov.SetRotation(hcovel);
1498 // now rotate in local system
1499 hcov.Multiply(svMatrix);
1500 hcov.MultiplyLeft(&svMatrix->Inverse());
6526a72c 1501 //
7b85e477 1502 // set local sigmas
1503 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
1504 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4])); // RS
1505 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
1506
1507 // set minimum value for SigmaLoc to 10 micron
1508 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
1509 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
6526a72c 1510 //
1511 /* RRR the rescaling is moved to PrepareTrack
7b85e477 1512 // multiply local sigmas by global and module specific factor
1513 for (int i=3;i--;) fSigmaLoc[i] *= fSigmaFactor[i]*fCurrentModule->GetSigmaFactor(i);
1514 //
6526a72c 1515 */
7b85e477 1516 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
1517
1518 return 0;
1519}
1520
6526a72c 1521//________________________________________________________________________________________________________
7b85e477 1522void AliITSAlignMille2::Print(Option_t*) const
1523{
45993b23 1524 // print current status
7b85e477 1525 printf("*** AliMillepede for ITS ***\n");
1526 printf(" Number of defined super modules: %d\n",fNModules);
1527 printf(" Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
1528 //
1529 if (fGeoManager)
1530 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
1531 else
1532 printf(" geometry not loaded\n");
1533 //
1534 if (fUsePreAlignment)
1535 printf(" using prealignment from %s \n",fPreAlignmentFileName.Data());
1536 else
1537 printf(" prealignment not used\n");
1538 //
1539 //
1540 if (fBOn)
1541 printf(" B Field set to %f T - using Riemann's helices\n",fBField);
1542 else
1543 printf(" B Field OFF - using straight lines \n");
1544 //
1545 if (fRequirePoints) printf(" Required points in tracks:\n");
1546 for (Int_t i=0; i<6; i++) {
1547 if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
1548 if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
1549 if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]);
1550 }
1551 for (Int_t i=0; i<3; i++) {
1552 if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
1553 if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
1554 if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]);
1555 }
1556 //
1557 printf("\n Millepede configuration parameters:\n");
7b85e477 1558 printf(" init value for chi2 cut : %.4f\n",fStartFac);
1559 printf(" first iteration cut value : %.4f\n",fResCutInitial);
1560 printf(" other iterations cut value : %.4f\n",fResCut);
1561 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
6526a72c 1562 printf(" def.scaling for local sigmas : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
7b85e477 1563
1564 printf("List of defined modules:\n");
1565 printf(" intidx\tindex\tvoluid\tname\n");
1566 for (int i=0; i<fNModules; i++) {
1567 AliITSAlignMille2Module* md = GetMilleModule(i);
1568 printf(" %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
1569 }
1570}
1571
6526a72c 1572//________________________________________________________________________________________________________
7b85e477 1573AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
1574{
6526a72c 1575 // return pointer to a defined supermodule
7b85e477 1576 // return NULL if error
6526a72c 1577 Int_t i=IsVIDDefined(voluid);
7b85e477 1578 if (i<0) return NULL;
1579 return GetMilleModule(i);
1580}
1581
6526a72c 1582//________________________________________________________________________________________________________
1583AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
7b85e477 1584{
6526a72c 1585 // return pointer to a defined supermodule
1586 // return NULL if error
1587 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1588 if (vid>0) return GetMilleModuleByVID(vid);
1589 else { // this is not alignable module, need to look within defined supermodules
1590 int i = IsSymDefined(symname);
1591 if (i>=0) return GetMilleModule(i);
1592 }
1593 return 0;
1594}
1595
1596//________________________________________________________________________________________________________
1597AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
1598{
1599 // return pointer to a defined/contained supermodule
1600 // return NULL otherwise
1601 int i = IsSymContained(symname);
1602 return i<0 ? 0 : GetMilleModule(i);
1603}
1604
1605//________________________________________________________________________________________________________
1606AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
1607{
45993b23 1608 // get delta from prealignment for given volume
6526a72c 1609 if (!fPrealignment) return 0;
1610 for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1611 AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
1612 if (!strcmp(preob->GetSymName(),symname)) return preob;
1613 }
1614 return 0;
1615}
1616
1617//________________________________________________________________________________________________________
1618AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
1619{
45993b23 1620 // get delta with respect to which the constraint is declared
6526a72c 1621 if (!fConstrRef) return 0;
1622 for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1623 AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
1624 if (!strcmp(preob->GetSymName(),symname)) return preob;
1625 }
1626 return 0;
1627}
1628
1629//________________________________________________________________________________________________________
1630Bool_t AliITSAlignMille2::InitRiemanFit()
1631{
1632 // Initialize Riemann Fitter for current track
1633 // return kFALSE if error
7b85e477 1634
1635 if (!fBOn) return kFALSE;
1636
1637 Int_t npts=0;
1638 AliTrackPoint ap;
1639 npts = fTrack->GetNPoints();
1640 AliDebug(3,Form("Fitting track with %d points",npts));
1641
1642 fRieman->Reset();
1643 fRieman->SetTrackPointArray(fTrack);
1644
1645 TArrayI ai(npts);
1646 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
1647
1648 // fit track with 5 params in his own tracking-rotated reference system
1649 // xc = -p[1]/p[0];
1650 // yc = 1/p[0];
1651 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
1652 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
1653 return kFALSE;
1654 }
1655
1656 for (int i=0; i<5; i++)
1657 fLocalInitParam[i] = fRieman->GetParam()[i];
1658
1659 return kTRUE;
1660}
1661
6526a72c 1662//________________________________________________________________________________________________________
45993b23 1663void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int flag)
7b85e477 1664{
45993b23 1665 // local function for minuit
7b85e477 1666 const double kTiny = 1.e-14;
1667 chi2 = 0;
1668 static AliTrackPoint pnt;
45993b23 1669 static Bool_t fullErr2D;
7b85e477 1670 //
45993b23 1671 if (flag==1) fullErr2D = kFALSE;//kTRUE;
7b85e477 1672 enum {kAX,kAZ,kBX,kBZ};
1673 enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1674 //
1675 AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
1676 AliTrackPointArray* track = alig->GetCurrentTrack();
1677 //
1678 int npts = track->GetNPoints();
1679 for (int ip=0;ip<npts;ip++) {
1680 track->GetPoint(pnt,ip);
1681 const float *cov = pnt.GetCov();
1682 double y = pnt.GetY();
1683 double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
1684 double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
1685 double xxe = cov[kXX];
1686 double zze = cov[kZZ];
1687 double xze = cov[kXZ];
1688 //
1689 if (fullErr2D) {
1690 xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
1691 zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
1692 xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
1693 }
1694 //
1695 double det = xxe*zze - xze*xze;
1696 if (det<kTiny) {
1697 printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
1698 "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
1699 xxe = cov[kXX];
1700 zze = cov[kZZ];
1701 xze = cov[kXZ];
1702 fullErr2D = kFALSE;
1703 }
1704 double xxeI = zze/det;
1705 double zzeI = xxe/det;
1706 double xzeI =-xze/det;
1707 //
1708 chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
1709 //
1710 // printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI, chi2);
1711 }
1712 //
1713}
1714
6526a72c 1715//________________________________________________________________________________________________________
7b85e477 1716void AliITSAlignMille2::InitTrackParams(int meth)
1717{
1718 /// initialize local parameters with different methods
1719 /// for current track (fTrack)
1720 Int_t npts=0;
1721 AliTrackPoint ap;
1722 double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
1723 // simple linear interpolation
1724 // get local starting parameters (to be substituted by ESD track parms)
1725 // local parms (fLocalInitParam[]) are:
1726 // [0] = global x coord. of straight line intersection at y=0 plane
1727 // [1] = global z coord. of straight line intersection at y=0 plane
1728 // [2] = px/py
1729 // [3] = pz/py
1730 // test #1: linear fit in x(y) and z(y)
1731 npts = fTrack->GetNPoints();
1732 AliDebug(3,Form("*** initializing track with %d points ***",npts));
1733 for (int i=npts;i--;) {
1734 sY += fTrack->GetY()[i];
1735 sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
1736 sX += fTrack->GetX()[i];
1737 sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
1738 sZ += fTrack->GetZ()[i];
1739 sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
1740 }
1741 det = sYY*npts-sY*sY;
1742 if (det==0) det = 1E-20;
1743 fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
1744 fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
1745 //
1746 fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
1747 fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
7b85e477 1748 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f ugx = %f\n",fLocalInitParam[0],fLocalInitParam[2]));
1749 //
1750 if (meth==1) return;
1751 //
1752 // perform full fit accounting for cov.matrix
1753 static TVirtualFitter *minuit = 0;
1754 static Double_t step[5] = {1E-3,1E-3,1E-4,1E-4,1E-5};
1755 static Double_t arglist[10];
1756 //
1757 if (!minuit) {
1758 minuit = TVirtualFitter::Fitter(0,4);
1759 minuit->SetFCN(trackFit2D);
1760 arglist[0] = 1;
1761 minuit->ExecuteCommand("SET ERR",arglist, 1);
1762 //
1763 arglist[0] = -1;
1764 minuit->ExecuteCommand("SET PRINT",arglist,1);
1765 //
1766 }
1767 //
1768 minuit->SetParameter(0, "ax", fLocalInitParam[0], step[0], 0,0);
1769 minuit->SetParameter(1, "az", fLocalInitParam[1], step[1], 0,0);
1770 minuit->SetParameter(2, "bx", fLocalInitParam[2], step[2], 0,0);
1771 minuit->SetParameter(3, "bz", fLocalInitParam[3], step[3], 0,0);
1772 //
1773 arglist[0] = 1000; // number of function calls
1774 arglist[1] = 0.001; // tolerance
7b85e477 1775 minuit->ExecuteCommand("MIGRAD",arglist,2);
1776 //
1777 for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
1778 for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
1779 //
1780}
1781
6526a72c 1782//________________________________________________________________________________________________________
1783Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
1784{
1785 // checks if supermodule with this symname is defined and return the internal index
1786 // return -1 if not.
1787 for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
1788 return -1;
1789}
1790
1791//________________________________________________________________________________________________________
1792Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
1793{
1794 // checks if module with this symname is defined and return the internal index
1795 // return -1 if not.
1796 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1797 if (vid>0) return IsVIDContained(vid);
1798 // only sensors have real vid, but maybe we have a supermodule with fake vid?
1799 // IMPORTANT: always start from the end to start from the sensors
1800 return IsSymDefined(symname);
1801}
7b85e477 1802
6526a72c 1803//________________________________________________________________________________________________________
1804Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
7b85e477 1805{
1806 // checks if supermodule 'voluid' is defined and return the internal index
6526a72c 1807 // return -1 if not.
1808 for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
1809 return -1;
7b85e477 1810}
1811
6526a72c 1812//________________________________________________________________________________________________________
1813Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
7b85e477 1814{
6526a72c 1815 // checks if the sensitive module 'voluid' is contained inside a supermodule
1816 // and return the internal index of the last identified supermodule
7b85e477 1817 // return -1 if error
6526a72c 1818 // IMPORTANT: always start from the end to start from the sensors
7b85e477 1819 if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
6526a72c 1820 for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
1821 return -1;
7b85e477 1822}
1823
6526a72c 1824//________________________________________________________________________________________________________
1825Int_t AliITSAlignMille2::CheckCurrentTrack()
7b85e477 1826{
7b85e477 1827 /// checks if AliTrackPoints belongs to defined modules
1828 /// return number of good poins
1829 /// return 0 if not enough points
1830
1831 Int_t npts = fTrack->GetNPoints();
1832 Int_t ngoodpts=0;
1833 // debug points
6526a72c 1834 for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
1835 //
7b85e477 1836 if (ngoodpts<fMinNPtsPerTrack) return 0;
1837
1838 return ngoodpts;
1839}
1840
6526a72c 1841//________________________________________________________________________________________________________
1842Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track)
1843{
7b85e477 1844 /// Process track; Loop over hits and set local equations
1845 /// here 'track' is a AliTrackPointArray
1846 /// return 0 if success;
1847
6526a72c 1848 if (!fIsMilleInit) Init();
1849 //
7b85e477 1850 Int_t npts = track->GetNPoints();
1851 AliDebug(2,Form("*** Input track with %d points ***",npts));
1852
1853 // preprocessing of the input track: keep only points in defined volumes,
1854 // move points if prealignment is set, sort by Yglo if required
1855
1856 fTrack=PrepareTrack(track);
1857 if (!fTrack) return -1;
1858
1859 npts = fTrack->GetNPoints();
6526a72c 1860 if (npts>kMaxPoints) {
1861 AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
1862 }
7b85e477 1863 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
1864
1865 if (!fBOn) { // straight lines
1866 // set local starting parameters (to be substituted by ESD track parms)
1867 // local parms (fLocalInitParam[]) are:
1868 // [0] = global x coord. of straight line intersection at y=0 plane
1869 // [1] = global z coord. of straight line intersection at y=0 plane
1870 // [2] = px/py
1871 // [3] = pz/py
1872 InitTrackParams(fInitTrackParamsMeth);
1873 }
1874 else {
1875 // local parms (fLocalInitParam[]) are the Riemann Fitter params
1876 if (!InitRiemanFit()) {
1877 AliInfo("Riemann fit failed! skipping this track...");
7b85e477 1878 fTrack=NULL;
1879 return -5;
1880 }
1881 }
1882
1883 Int_t nloceq=0;
6526a72c 1884 Int_t ngloeq=0;
1885 static Mille2Data md[kMaxPoints];
1886 //
7b85e477 1887 for (Int_t ipt=0; ipt<npts; ipt++) {
1888 fTrack->GetPoint(fCluster,ipt);
6526a72c 1889 fCluster.SetUniqueID(ipt);
7b85e477 1890 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
1891
1892 // set geometry parameters for the the current module
1893 if (InitModuleParams()) continue;
1894 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n",
1895 track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
1896 fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
1897 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
6526a72c 1898 int res = AddLocalEquation(md[nloceq]);
1899 if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
1900 else if (res==0) nloceq++;
1901 else {nloceq++; ngloeq++;}
7b85e477 1902 } // end loop over points
1903 //
7b85e477 1904 fTrack=NULL;
6526a72c 1905 // not enough good points?
1906 if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
7b85e477 1907 //
1908 // finally send local equations to millepede
1909 SetLocalEquations(md,nloceq);
1910 fMillepede->SaveRecordData(); // RRR
1911 //
1912 return 0;
1913}
1914
6526a72c 1915//________________________________________________________________________________________________________
1916Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar)
1917{
7b85e477 1918 /// calculate track intersection point in local coordinates
1919 /// according with a given set of parameters (local(4) and global(6))
1920 /// and fill fPintLoc/Glo
1921 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
1922 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
1923 /// return 0 if success
1924
1925 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]));
1926 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]));
1927
1928
1929 // prepare the TGeoHMatrix
6526a72c 1930 TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
1931 !fUseGlobalDelta);
1932 if (!tempHMat) return -1;
7b85e477 1933
1934 Double_t v0g[3]; // vector with straight line direction in global coord.
1935 Double_t p0g[3]; // point of the straight line (glo)
1936
1937 if (fBOn) { // B FIELD!
1938 AliTrackPoint prf;
1939 for (int ip=0; ip<5; ip++)
1940 fRieman->SetParam(ip,lpar[ip]);
1941
1942 if (!fRieman->GetPCA(fCluster,prf)) {
1943 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
1944 return -3;
1945 }
1946 // now determine straight line passing tangent to fit curve at prf
1947 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
1948 // mo' P1=(X,Y,Z)_glo_prf
1949 // => (x,y,Z)_trk_prf ruotando di alpha...
1950 Double_t alpha=fRieman->GetAlpha();
1951 Double_t x1g = prf.GetX();
1952 Double_t y1g = prf.GetY();
1953 Double_t z1g = prf.GetZ();
1954 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
1955 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
1956 Double_t z1t = z1g;
1957
1958 Double_t x2t = x1t+1.0;
1959 Double_t y2t = y1t+fRieman->GetDYat(x1t);
1960 Double_t z2t = z1t+fRieman->GetDZat(x1t);
1961 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
1962 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
1963 Double_t z2g = z2t;
1964
1965 AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi()));
1966 AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
1967 AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
1968
1969 if (TMath::Abs(y2g-y1g)<1e-15) {
1970 AliInfo("DeltaY=0! Cannot proceed...");
1971 return -1;
1972 }
1973 // ugx,1,ugz
1974 v0g[0] = (x2g-x1g)/(y2g-y1g);
1975 v0g[1]=1.0;
1976 v0g[2] = (z2g-z1g)/(y2g-y1g);
1977
1978 // point: just keep prf
1979 p0g[0]=x1g;
1980 p0g[1]=y1g;
1981 p0g[2]=z1g;
1982 }
1983 else { // staight line
1984 // vector of initial straight line direction in glob. coord
1985 v0g[0]=lpar[2];
1986 v0g[1]=1.0;
1987 v0g[2]=lpar[3];
1988
1989 // intercept in yg=0 plane in glob coord
1990 p0g[0]=lpar[0];
1991 p0g[1]=0.0;
1992 p0g[2]=lpar[1];
1993 }
1994 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]));
1995
1996 // same in local coord.
1997 Double_t p0l[3],v0l[3];
6526a72c 1998 tempHMat->MasterToLocalVect(v0g,v0l);
1999 tempHMat->MasterToLocal(p0g,p0l);
7b85e477 2000
2001 if (TMath::Abs(v0l[1])<1e-15) {
2002 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2003 return -1;
2004 }
2005
2006 // local intersection point
2007 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2008 fPintLoc[1] = 0;
2009 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2010
2011 // global intersection point
6526a72c 2012 tempHMat->LocalToMaster(fPintLoc,fPintGlo);
7b85e477 2013 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]));
2014
2015 return 0;
2016}
2017
6526a72c 2018//________________________________________________________________________________________________________
2019Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
2020{
7b85e477 2021 /// calculate numerically (ROOT's style) the derivatives for
2022 /// local X intersection and local Z intersection
2023 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
2024 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2025 /// return 0 if success
2026
2027 // copy initial parameters
6526a72c 2028 Double_t lpar[kNLocal];
2029 Double_t gpar[kNParCh];
2030 Double_t *derivative;
2031 for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2032 for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
7b85e477 2033
2034 // trial with fixed dpar...
6526a72c 2035 Double_t dpar = 0.;
7b85e477 2036
2037 if (islpar) { // track parameters
2038 //dpar=fLocalInitParam[paridx]*0.001;
2039 // set minimum dpar
6526a72c 2040 derivative = fDerivativeLoc[paridx];
7b85e477 2041 if (!fBOn) {
6526a72c 2042 if (paridx<3) dpar=1.0e-4; // translations
7b85e477 2043 else dpar=1.0e-6; // direction
2044 }
2045 else { // B Field
2046 // pepo: proviamo con 1/1000, poi evenctually 1/100...
2047 Double_t dfrac=0.01;
2048 switch(paridx) {
2049 case 0:
2050 // RMS cosmics: 1e-4
2051 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2052 break;
2053 case 1:
2054 // RMS cosmics: 0.2
2055 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2056 break;
2057 case 2:
2058 // RMS cosmics: 9
2059 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2060 break;
2061 case 3:
2062 // RMS cosmics: 7
2063 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2064 break;
2065 case 4:
2066 // RMS cosmics: 0.3
2067 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2068 break;
2069 }
2070 }
2071 }
2072 else { // alignment global parameters
6526a72c 2073 derivative = fDerivativeGlo[paridx];
7b85e477 2074 //dpar=fModuleInitParam[paridx]*0.001;
2075 if (paridx<3) dpar=1.0e-4; // translations
2076 else dpar=1.0e-2; // angles
2077 }
2078
2079 AliDebug(3,Form("+++ using dpar=%g",dpar));
2080
2081 // calculate derivative ROOT's like:
2082 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2083 Double_t pintl1[3]; // f(x-h)
2084 Double_t pintl2[3]; // f(x-h/2)
2085 Double_t pintl3[3]; // f(x+h/2)
2086 Double_t pintl4[3]; // f(x+h)
2087
2088 // first values
2089 if (islpar) lpar[paridx] -= dpar;
2090 else gpar[paridx] -= dpar;
2091 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2092 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2093
2094 // second values
2095 if (islpar) lpar[paridx] += dpar/2;
2096 else gpar[paridx] += dpar/2;
2097 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2098 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2099
2100 // third values
2101 if (islpar) lpar[paridx] += dpar;
2102 else gpar[paridx] += dpar;
2103 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2104 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2105
2106 // fourth values
2107 if (islpar) lpar[paridx] += dpar/2;
2108 else gpar[paridx] += dpar/2;
2109 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2110 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2111
2112 Double_t h2 = 1./(2.*dpar);
2113 Double_t d0 = pintl4[0]-pintl1[0];
2114 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
6526a72c 2115 derivative[0] = h2*(4*d2 - d0)/3.;
2116 if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
7b85e477 2117
2118 d0 = pintl4[2]-pintl1[2];
2119 d2 = 2.*(pintl3[2]-pintl2[2]);
6526a72c 2120 derivative[2] = h2*(4*d2 - d0)/3.;
2121 if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
7b85e477 2122
2123 AliDebug(3,Form("\n+++ derivatives +++ \n"));
6526a72c 2124 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2125 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
7b85e477 2126
2127 return 0;
2128}
2129
6526a72c 2130//________________________________________________________________________________________________________
2131Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2132{
7b85e477 2133 /// Define local equation for current cluster in X and Z coor.
2134 /// and store them to memory
6526a72c 2135 /// return -1 in case of failure to build some equation
2136 /// 0 if no free global parameters were found but local eq is built
2137 /// 1 if both local and global eqs are built
2138 //
45993b23 2139 // store first intersection point
6526a72c 2140 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
7b85e477 2141 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2142 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2143
2144 // calculate local derivatives numerically
2145 Bool_t zeroX = kTRUE;
2146 Bool_t zeroZ = kTRUE;
2147 //
2148 for (Int_t i=0; i<fNLocal; i++) {
2149 if (CalcDerivatives(i,kTRUE)) return -1;
45993b23 2150 m.fDerLocX[i] = fDerivativeLoc[i][0];
2151 m.fDerLocZ[i] = fDerivativeLoc[i][2];
6526a72c 2152 if (zeroX) zeroX = fDerivativeLoc[i][0]==0;
2153 if (zeroZ) zeroZ = fDerivativeLoc[i][2]==0;
7b85e477 2154 }
2155 // for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
7b85e477 2156 //
45993b23 2157 if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2158 if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
6526a72c 2159 //
2160 int ifill = 0;
7b85e477 2161 //
2162 AliITSAlignMille2Module* endModule = fCurrentModule;
2163 //
6526a72c 2164 zeroX = zeroZ = kTRUE;
2165 Bool_t dfDone[kNParCh];
2166 for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
45993b23 2167 m.fNModFilled = 0;
6526a72c 2168 //
2169 // special block for SDD derivatives
2170 Double_t jacobian[kNParChGeom];
2171 Int_t nmodTested = 0;
2172 //
7b85e477 2173 do {
6526a72c 2174 if (fCurrentModule->GetNParFree()==0) continue;
2175 nmodTested++;
2176 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2177 //
2178 if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2179 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2180 if (!dfDone[i]) {
2181 if (CalcDerivatives(i,kFALSE)) return -1;
2182 else {
2183 dfDone[i] = kTRUE;
2184 if (zeroX) zeroX = fDerivativeGlo[i][0]==0;
2185 if (zeroZ) zeroZ = fDerivativeGlo[i][2]==0;
2186 }
2187 }
2188 //
45993b23 2189 m.fDerGloX[ifill] = fDerivativeGlo[i][0];
2190 m.fDerGloZ[ifill] = fDerivativeGlo[i][2];
2191 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
6526a72c 2192 }
2193 //
2194 // specific for special sensors
2195 if ( fCurrentModule->IsSDD() &&
2196 (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0 ||
2197 fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) ) {
2198 //
2199 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2200 // where V0 and T are the nominal drift velocity, time and time0
2201 // and the dT0 and dV are the corrections:
2202 // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2203 // dX/dV = dX/dxloc * dxloc/dV = dX/dxloc * (T-T0)
2204 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2205 //
2206 if (!dfDone[AliITSAlignMille2Module::kDOFT0] || !dfDone[AliITSAlignMille2Module::kDOFDV]) {
2207 //
2208 double dXdxlocsens=0., dZdxlocsens=0.;
2209 //
2210 // if the current module is the sensor itself and we work with local params, then
2211 // we can directly take dX/dxloc_sens dZ/dxloc_sens
2212 if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2213 if (dfDone[AliITSAlignMille2Module::kDOFTX]) {
2214 CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2215 dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2216 }
2217 dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2218 dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2219 }
2220 //
2221 else { // need to perform some transformations
2222 // fetch the jacobian of the transformation from the sensors local frame to the frame
2223 // where the parameters are defined:
2224 // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2225 if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2226 AliITSAlignMille2Module::kDOFTX, jacobian);
2227 // Local: dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2228 else fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2229 AliITSAlignMille2Module::kDOFTX, jacobian);
2230 //
2231 for (int j=0;j<kNParChGeom;j++) {
2232 // need global derivative even if the j-th param is locked
2233 if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2234 dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2235 dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2236 }
2237 }
2238 //
2239 if (zeroX) zeroX = dXdxlocsens == 0;
2240 if (zeroZ) zeroZ = dZdxlocsens == 0;
2241 //
2242 double vdrift = GetVDriftSDD();
2243 double tdrift = GetTDriftSDD();
2244 //
2245 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2246 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2247 dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2248 //
2249 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0] = -dXdxlocsens*TMath::Sign(tdrift,vdrift);
2250 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2] = -dZdxlocsens*TMath::Sign(tdrift,vdrift);
2251 dfDone[AliITSAlignMille2Module::kDOFDV] = kTRUE;
2252 //
2253 }
2254 //
2255 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
45993b23 2256 m.fDerGloX[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2257 m.fDerGloZ[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2258 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
6526a72c 2259 }
2260 //
2261 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) {
45993b23 2262 m.fDerGloX[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0];
2263 m.fDerGloZ[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2];
2264 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV);
7b85e477 2265 }
2266 }
7b85e477 2267 //
45993b23 2268 m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
7b85e477 2269 } while( (fCurrentModule=fCurrentModule->GetParent()) );
2270 //
45993b23 2271 if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2272 if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
6526a72c 2273 //
7b85e477 2274 // ok, can copy to m
2275 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
45993b23 2276 m.fMeasX = fMeasLoc[0]-fPintLoc0[0];
2277 m.fSigmaX = fSigmaLoc[0];
7b85e477 2278 //
2279 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
45993b23 2280 m.fMeasZ = fMeasLoc[2]-fPintLoc0[2];
2281 m.fSigmaZ = fSigmaLoc[2];
7b85e477 2282 //
45993b23 2283 m.fNGlobFilled = ifill;
7b85e477 2284 fCurrentModule = endModule;
2285 //
6526a72c 2286 return Int_t(!zeroX && !zeroZ);
7b85e477 2287}
2288
6526a72c 2289//________________________________________________________________________________________________________
2290void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
2291{
7b85e477 2292 /// Set local equations with data stored in m
2293 /// return 0 if success
2294 //
2295 for (Int_t j=0; j<neq; j++) {
2296 //
2297 const Mille2Data &m = marr[j];
2298 //
2299 // set equation for Xloc coordinate
45993b23 2300 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",m.fMeasX, m.fSigmaX));
2301 for (int i=fNLocal; i--;) SetLocalDerivative( i, m.fDerLocX[i] );
2302 for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGloX[i] );
2303 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeasX, m.fSigmaX);
7b85e477 2304 //
2305 // set equation for Zloc coordinate
45993b23 2306 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",m.fMeasZ, m.fSigmaZ));
2307 for (int i=fNLocal; i--;) SetLocalDerivative( i, m.fDerLocZ[i] );
2308 for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGloZ[i] );
2309 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeasZ, m.fSigmaZ);
6526a72c 2310 //
45993b23 2311 for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
6526a72c 2312 //
7b85e477 2313 }
2314}
2315
6526a72c 2316//________________________________________________________________________________________________________
2317Int_t AliITSAlignMille2::GlobalFit()
2318{
7b85e477 2319 /// Call global fit; Global parameters are stored in parameters
6526a72c 2320 if (!fIsMilleInit) Init();
2321 //
2322 ApplyPreConstraints();
2323 int res = fMillepede->GlobalFit();
7b85e477 2324 AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
6526a72c 2325 if (res) {
2326 // fetch the parameters
2327 for (int imd=fNModules;imd--;) {
2328 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2329 int nprocp = 0;
2330 for (int ip=mod->GetNParTot();ip--;) {
2331 int idp = mod->GetParOffset(ip);
2332 if (idp<0) continue; // was not in the explicit fit
2333 mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
2334 mod->SetParErr(ip,fMillepede->GetFinalError(idp));
2335 int np = fMillepede->GetProcessedPoints(idp);
2336 if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
2337 }
2338 if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
2339 }
7b85e477 2340
7b85e477 2341 }
6526a72c 2342 ApplyPostConstraints();
2343 return res;
7b85e477 2344}
2345
6526a72c 2346//________________________________________________________________________________________________________
2347void AliITSAlignMille2::PrintGlobalParameters()
2348{
7b85e477 2349 /// Print global parameters
2350 if (!fIsMilleInit) {
2351 AliInfo("Millepede has not been initialized!");
2352 return;
2353 }
2354 fMillepede->PrintGlobalParameters();
2355}
2356
6526a72c 2357//________________________________________________________________________________________________________
7b85e477 2358Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
2359{
2360 // load definitions of supermodules from a root file
2361 // return 0 if success
2362
2363 TFile *smf=TFile::Open(sfile);
2364 if (!smf->IsOpen()) {
2365 AliInfo(Form("Cannot open supermodule file %s",sfile));
2366 return -1;
2367 }
2368
2369 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
2370 if (!sma) {
2371 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
2372 return -2;
2373 }
2374 Int_t nsma=sma->GetEntriesFast();
2375 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
2376 //
2377 Char_t st[250];
2378 char symname[150];
2379 UShort_t volid;
2380 TGeoHMatrix m;
2381 //
2382 for (Int_t i=0; i<nsma; i++) {
2383 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
2384 volid=a->GetVolUID();
2385 strcpy(st,a->GetSymName());
2386 a->GetMatrix(m);
2387 //
2388 sscanf(st,"%s",symname);
6526a72c 2389 //
7b85e477 2390 // decode module list
2391 char *stp=strstr(st,"ModuleList:");
2392 if (!stp) return -3;
2393 stp += 11;
2394 int idx[2200];
2395 char spp[200]; int jp=0;
2396 char cl[20];
2397 strcpy(st,stp);
2398 int l=strlen(st);
2399 int j=0;
2400 int n=0;
2401 //
2402 while (j<=l) {
2403 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
2404 spp[jp]=0;
2405 jp=0;
2406 if (strlen(spp)) {
2407 int k=strcspn(spp,"-");
2408 if (k<int(strlen(spp))) { // c'e' il -
2409 strcpy(cl,&(spp[k+1]));
2410 spp[k]=0;
2411 int ifrom=atoi(spp); int ito=atoi(cl);
2412 for (int b=ifrom; b<=ito; b++) {
2413 idx[n]=b;
2414 n++;
2415 }
2416 }
2417 else { // numerillo singolo
2418 idx[n]=atoi(spp);
2419 n++;
2420 }
2421 }
2422 }
2423 else {
2424 spp[jp]=st[j];
2425 jp++;
2426 }
2427 j++;
2428 }
2429 UShort_t volidsv[2198];
2430 for (j=0;j<n;j++) {
2431 volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
2432 if (!volidsv[j]) {
2433 AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
2434 return -5;
2435 }
2436 }
2437 Int_t smindex=int(2198+volid-14336); // virtual index
6526a72c 2438 //
7b85e477 2439 fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
2440 //
2441 fNSuperModules++;
2442 }
6526a72c 2443 //
7b85e477 2444 smf->Close();
2445 //
2446 return 0;
2447}
2448
6526a72c 2449//________________________________________________________________________________________________________
2450void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2451{
2452 // require that sum of modifications for the childs of this module is = val, i.e.
2453 // the internal corrections moves the module as a whole by fixed value (0 by default).
2454 // pattern is the bit pattern for the parameters to constrain
2455 //
6526a72c 2456 if (fIsMilleInit) {
2457 AliInfo("Millepede has been already initialized: no constrain may be added!");
2458 return;
7b85e477 2459 }
6526a72c 2460 if (!GetMilleModule(idm)->GetNChildren()) return;
2461 TString nm = "cstrSUMean";
2462 nm += GetNConstraints();
2463 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2464 idm,val,pattern);
2465 cstr->SetConstraintID(GetNConstraints());
2466 fConstraints.Add(cstr);
7b85e477 2467}
2468
6526a72c 2469//________________________________________________________________________________________________________
2470void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2471{
2472 // require that median of the modifications for the childs of this module is = val, i.e.
2473 // the internal corrections moves the module as a whole by fixed value (0 by default)
2474 // module the outliers.
2475 // pattern is the bit pattern for the parameters to constrain
2476 // The difference between the mean and the median will be transfered to the parent
6526a72c 2477 if (fIsMilleInit) {
2478 AliInfo("Millepede has been already initialized: no constrain may be added!");
2479 return;
2480 }
2481 if (!GetMilleModule(idm)->GetNChildren()) return;
2482 TString nm = "cstrSUMed";
2483 nm += GetNConstraints();
2484 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2485 idm,val,pattern);
2486 cstr->SetConstraintID(GetNConstraints());
2487 fConstraints.Add(cstr);
2488}
2489
2490//________________________________________________________________________________________________________
2491void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
2492{
2493 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2494 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2495 // pattern is the bit pattern for the parameters to constrain
7b85e477 2496 //
6526a72c 2497 if (fIsMilleInit) {
2498 AliInfo("Millepede has been already initialized: no constrain may be added!");
2499 return;
2500 }
2501 TString nm = "cstrOMean";
2502 nm += GetNConstraints();
2503 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2504 -1,val,pattern);
2505 cstr->SetConstraintID(GetNConstraints());
2506 fConstraints.Add(cstr);
2507}
2508
2509//________________________________________________________________________________________________________
2510void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
2511{
2512 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2513 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2514 // pattern is the bit pattern for the parameters to constrain
7b85e477 2515 //
6526a72c 2516 if (fIsMilleInit) {
2517 AliInfo("Millepede has been already initialized: no constrain may be added!");
2518 return;
7b85e477 2519 }
6526a72c 2520 TString nm = "cstrOMed";
2521 nm += GetNConstraints();
2522 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2523 -1,val,pattern);
2524 cstr->SetConstraintID(GetNConstraints());
2525 fConstraints.Add(cstr);
2526}
2527
2528//________________________________________________________________________________________________________
2529void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
2530{
45993b23 2531 // apply constraint on parameters in the local frame
6526a72c 2532 if (fIsMilleInit) {
2533 AliInfo("Millepede has been already initialized: no constrain may be added!");
2534 return;
2535 }
2536 AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
2537 cstr->SetConstraintID(GetNConstraints());
2538 fConstraints.Add(cstr);
2539}
2540
2541//________________________________________________________________________________________________________
45993b23 2542void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
6526a72c 2543{
2544 // apply the constraint on the local corrections of a list of modules
2545 int nmod = cstr->GetNModules();
2546 double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
7b85e477 2547 //
6526a72c 2548 for (int imd=nmod;imd--;) {
2549 int modID = cstr->GetModuleID(imd);
2550 AliITSAlignMille2Module* mod = GetMilleModule(modID);
2551 ResetLocalEquation();
2552 int nadded = 0;
2553 double value = cstr->GetValue();
2554 double sigma = cstr->GetError();
2555 //
2556 // in case the reference (survey) deltas were imposed for Gaussian constraints
2557 // already accumulated corrections: they must be subtracted from the constraint value.
2558 if (IsConstraintWrtRef()) {
2559 //
2560 Double_t precal[AliITSAlignMille2Module::kMaxParTot];
2561 Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
2562 for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
2563 //
2564 // check if there was a reference delta provided for this module
2565 AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
2566 if (parref) parref->GetPars(refcal, refcal+3); // found reference delta
2567 //
2568 // extract already applied local corrections for this module
2569 if (fPrealignment) {
2570 //
2571 AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
2572 if (preo) {
2573 TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); // Delta_Glob * Delta_Glob_Par * M
2574 preo->GetMatrix(preMat); // Delta_Glob
2575 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
2576 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
2577 AliAlignObjParams algob;
2578 algob.SetMatrix(tmpMat);
2579 algob.GetPars(precal,precal+3); // local corrections for geometry
2580 }
2581 }
2582 //
2583 // subtract the contribution to constraint from precalibration
2584 for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
2585 //
2586 }
2587 //
2588 if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
2589 //
2590 for (int ipar=cstr->GetNCoeffs();ipar--;) {
2591 double coef = cstr->GetCoeff(ipar);
2592 if (coef==0) continue;
2593 //
2594 if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
2595 // we are working with local params or if the given param is not related to geometry,
2596 // apply the constraint directly
2597 int parPos = mod->GetParOffset(ipar);
2598 if (parPos<0) continue; // not in the fit
2599 fGlobalDerivatives[parPos] += coef;
2600 nadded++;
2601 }
2602 else { // we are working with global params, while the constraint is on local ones -> jacobian
2603 for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
2604 int parPos = mod->GetParOffset(jpar);
2605 if (parPos<0) continue;
2606 fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
2607 nadded++;
2608 }
2609 }
2610 }
2611 if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
2612 }
7b85e477 2613 //
6526a72c 2614}
2615
2616//________________________________________________________________________________________________________
2617void AliITSAlignMille2::ApplyPreConstraints()
2618{
45993b23 2619 // apply constriants which cannot be imposed after the fit
6526a72c 2620 int nconstr = GetNConstraints();
2621 for (int i=0;i<nconstr;i++) {
2622 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
7b85e477 2623 //
6526a72c 2624 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
2625 ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
2626 continue;
2627 }
7b85e477 2628 //
6526a72c 2629 if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
2630 //
2631 if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
2632 // apply constraint on the mean's before the fit
2633 int imd = cstr->GetModuleID();
2634 if (imd>=0) {
2635 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2636 UInt_t pattern = 0;
2637 for (int ipar=mod->GetNParTot();ipar--;) {
2638 if (!cstr->IncludesParam(ipar)) continue;
2639 if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
2640 pattern |= 0x1<<ipar;
2641 cstr->SetApplied(ipar);
2642 }
2643 ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
2644 //
2645 }
2646 else if (!PseudoParentsAllowed()) {
2647 ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
2648 cstr->SetApplied(-1);
2649 }
7b85e477 2650 }
6526a72c 2651}
2652
2653//________________________________________________________________________________________________________
2654void AliITSAlignMille2::ApplyPostConstraints()
2655{
45993b23 2656 // apply constraints which can be imposed after the fit
6526a72c 2657 int nconstr = GetNConstraints();
2658 Bool_t convGlo = kFALSE;
2659 // check if there is something to do
2660 int ntodo = 0;
2661 for (int i=0;i<nconstr;i++) {
2662 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2663 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2664 if (cstr->GetRemainingPattern() == 0) continue;
2665 ntodo++;
2666 }
2667 if (!ntodo) return;
7b85e477 2668 //
6526a72c 2669 if (!fUseGlobalDelta) { // need to convert to global params
2670 ConvertParamsToGlobal();
2671 convGlo = kTRUE;
2672 }
7b85e477 2673 //
6526a72c 2674 for (int i=0;i<nconstr;i++) {
2675 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2676 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2677 //
2678 int imd = cstr->GetModuleID();
2679 //
2680 if (imd>=0) {
2681 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2682 UInt_t pattern = 0;
2683 for (int ipar=mod->GetNParTot();ipar--;) {
2684 if (cstr->IsApplied(ipar)) continue;
2685 if (!cstr->IncludesParam(ipar)) continue;
2686 if (!mod->IsFreeDOF(ipar)) continue; // parameter is fixed, will not apply constraint
2687 pattern |= 0x1<<ipar;
2688 cstr->SetApplied(ipar);
2689 }
2690 if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
2691 //
2692 }
2693 else if (PseudoParentsAllowed()) {
2694 UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
2695 PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
2696 cstr->SetApplied(-1);
2697 }
2698 }
2699 // if there was a conversion, rewind it
2700 if (convGlo) ConvertParamsToLocal();
2701 //
7b85e477 2702}
2703
6526a72c 2704//________________________________________________________________________________________________________
2705void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2706{
6526a72c 2707 // require that sum of modifications for the childs of this module is = val, i.e.
2708 // the internal corrections moves the module as a whole by fixed value (0 by default).
7b85e477 2709 // pattern is the bit pattern for the parameters to constrain
2710 //
7b85e477 2711 //
6526a72c 2712 AliITSAlignMille2Module* mod = GetMilleModule(idm);
7b85e477 2713 //
6526a72c 2714 for (int ip=0;ip<kNParCh;ip++) {
2715 if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
2716 ResetLocalEquation();
2717 int nadd = 0;
2718 for (int ich=mod->GetNChildren();ich--;) {
2719 int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
2720 if (idpar<0) continue;
2721 fGlobalDerivatives[idpar] = 1.0;
2722 nadd++;
2723 }
2724 //
2725 if (nadd>0) {
2726 AddConstraint(fGlobalDerivatives,val);
2727 AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
2728 }
7b85e477 2729 }
7b85e477 2730 //
6526a72c 2731}
2732
2733//________________________________________________________________________________________________________
2734void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
2735{
2736 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2737 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2738 // pattern is the bit pattern for the parameters to constrain
2739 //
7b85e477 2740 for (int ip=0;ip<kNParCh;ip++) {
6526a72c 2741 //
2742 if ( !((pattern>>ip)&0x1) ) continue;
7b85e477 2743 ResetLocalEquation();
6526a72c 2744 int nadd = 0;
2745 for (int imd=fNModules;imd--;) {
2746 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2747 if (mod->GetParent()) continue; // this is not an orphan
2748 int idpar = mod->GetParOffset(ip);
2749 if (idpar<0) continue;
2750 fGlobalDerivatives[idpar] = 1.0;
2751 nadd++;
2752 }
2753 if (nadd>0) {
2754 AddConstraint(fGlobalDerivatives,val);
2755 AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
7b85e477 2756 }
7b85e477 2757 }
2758 //
7b85e477 2759 //
2760}
6526a72c 2761
2762//________________________________________________________________________________________________________
2763void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2764{
6526a72c 2765 // require that median or mean of the modifications for the childs of this module is = val, i.e.
2766 // the internal corrections moves the module as a whole by fixed value (0 by default)
2767 // module the outliers.
7b85e477 2768 // pattern is the bit pattern for the parameters to constrain
6526a72c 2769 // The difference between the mean and the median will be transfered to the parent
7b85e477 2770 //
6526a72c 2771 AliITSAlignMille2Module* parent = GetMilleModule(idm);
2772 int nc = parent->GetNChildren();
7b85e477 2773 //
6526a72c 2774 double *tmpArr = new double[nc];
7b85e477 2775 //
2776 for (int ip=0;ip<kNParCh;ip++) {
6526a72c 2777 int npc = 0;
2778 if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
2779 // compute the mean and median of the deltas
2780 int nfree = 0;
2781 for (int ich=nc;ich--;) {
2782 AliITSAlignMille2Module* child = parent->GetChild(ich);
2783 // if (!child->IsFreeDOF(ip)) continue;
2784 tmpArr[nfree++] = child->GetParVal(ip);
2785 }
2786 double median=0,mean=0;
2787 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
2788 mean += tmpArr[ic0];
2789 for (int ic1=ic0+1;ic1<nfree;ic1++)
2790 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
2791 }
7b85e477 2792 //
6526a72c 2793 int kmed = nfree/2;
2794 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
2795 if (nfree>0) mean /= nfree;
2796 //
2797 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
2798 //
2799 for (int ich=nc;ich--;) {
2800 AliITSAlignMille2Module* child = parent->GetChild(ich);
2801 // if (!child->IsFreeDOF(ip)) continue;
2802 child->SetParVal(ip, child->GetParVal(ip) + shift);
2803 npc++;
7b85e477 2804 }
6526a72c 2805 //
2806 parent->SetParVal(ip, parent->GetParVal(ip) - shift);
2807 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d children of module %d: %s",
2808 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
2809 ip,npc,idm,parent->GetName()));
7b85e477 2810 }
6526a72c 2811 delete[] tmpArr;
7b85e477 2812 //
7b85e477 2813 //
2814}
2815
6526a72c 2816//________________________________________________________________________________________________________
2817void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
7b85e477 2818{
6526a72c 2819 // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
2820 // the corrections moves the whole setup by fixed value (0 by default).
2821 // pattern is the bit pattern for the parameters to constrain
7b85e477 2822 //
6526a72c 2823 int nc = fNModules;
7b85e477 2824 //
6526a72c 2825 int norph = 0;
2826 for (int ich=nc;ich--;) if (!GetMilleModule(ich)->GetParent()) norph ++;
2827 if (!norph) return;
2828 double *tmpArr = new double[norph];
7b85e477 2829 //
6526a72c 2830 for (int ip=0;ip<kNParCh;ip++) {
2831 int npc = 0;
2832 if ( !((pattern>>ip)&0x1)) continue;
2833 // compute the mean and median of the deltas
2834 int nfree = 0;
2835 for (int ich=nc;ich--;) {
2836 AliITSAlignMille2Module* child = GetMilleModule(ich);
2837 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
2838 if (child->GetParent()) continue;
2839 tmpArr[nfree++] = child->GetParVal(ip);
2840 }
2841 double median=0,mean=0;
2842 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
2843 mean += tmpArr[ic0];
2844 for (int ic1=ic0+1;ic1<nfree;ic1++)
2845 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
2846 }
2847 //
2848 int kmed = nfree/2;
2849 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
2850 if (nfree>0) mean /= nfree;
2851 //
2852 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
2853 //
2854 for (int ich=nc;ich--;) {
2855 AliITSAlignMille2Module* child = GetMilleModule(ich);
2856 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
2857 if (child->GetParent()) continue;
2858 child->SetParVal(ip, child->GetParVal(ip) + shift);
2859 npc++;
7b85e477 2860 }
6526a72c 2861 //
2862 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d orphan modules",
2863 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
2864 ip,npc));
7b85e477 2865 }
6526a72c 2866 delete[] tmpArr;
7b85e477 2867 //
6526a72c 2868}
2869
2870//________________________________________________________________________________________________________
45993b23 2871Bool_t AliITSAlignMille2::IsParModConstrained(const AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
6526a72c 2872{
2873 // check if par of the module participates in some constraint, and set the flag for their types
2874 meanmed = gaussian = kFALSE;
2875 //
2876 if ( mod->IsParConstrained(par) ) gaussian = kTRUE; // direct constraint on this param
2877 //
2878 for (int icstr=GetNConstraints();icstr--;) {
2879 AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
2880 //
2881 if (!cstr->IncludesModPar(mod,par)) continue;
2882 if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
2883 else meanmed = kTRUE;
2884 //
2885 if (meanmed && gaussian) break; // no sense to check further
7b85e477 2886 }
6526a72c 2887 //
2888 return meanmed||gaussian;
2889}
2890
2891//________________________________________________________________________________________________________
45993b23 2892Bool_t AliITSAlignMille2::IsParModFamilyVaried(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
6526a72c 2893{
2894 // check if parameter par is varied for this module or its children up to the level depth
2895 if (depth<0) return kFALSE;
2896 if (mod->GetParOffset(par)>=0) return kTRUE;
2897 for (int icld=mod->GetNChildren();icld--;) {
2898 AliITSAlignMille2Module* child = mod->GetChild(icld);
2899 if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
7b85e477 2900 }
6526a72c 2901 return kFALSE;
7b85e477 2902 //
6526a72c 2903}
2904
2905/*
2906//________________________________________________________________________________________________________
2907Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
2908{
2909 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
2910 if (depth<0) return kTRUE;
2911 for (int icld=mod->GetNChildren();icld--;) {
2912 AliITSAlignMille2Module* child = mod->GetChild(icld);
2913 //if (child->GetParOffset(par)<0) continue; // fixed
2914 Bool_t cstMM=kFALSE,cstGS=kFALSE;
2915 // does this child have gaussian constraint ?
2916 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
2917 // check its children
2918 if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
2919 }
2920 return kFALSE;
2921 //
2922}
2923*/
2924
2925//________________________________________________________________________________________________________
45993b23 2926Bool_t AliITSAlignMille2::IsParFamilyFree(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
6526a72c 2927{
2928 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
2929 if (depth<0) return kFALSE;
2930 for (int icld=mod->GetNChildren();icld--;) {
2931 AliITSAlignMille2Module* child = mod->GetChild(icld);
2932 //if (child->GetParOffset(par)<0) continue; // fixed
2933 Bool_t cstMM=kFALSE,cstGS=kFALSE;
2934 // does this child have gaussian constraint ?
2935 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
2936 // check its children
2937 if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
2938 }
2939 return kFALSE;
7b85e477 2940 //
2941}
6526a72c 2942
2943//________________________________________________________________________________________________________
2944Double_t AliITSAlignMille2::GetTDriftSDD() const
2945{
45993b23 2946 // obtain drift time corrected for t0
6526a72c 2947 double t = fCluster.GetDriftTime();
2948 return t - fDriftTime0[ fCluster.GetUniqueID() ];
2949}
2950
2951//________________________________________________________________________________________________________
2952Double_t AliITSAlignMille2::GetVDriftSDD() const
2953{
45993b23 2954 // obtain corrected drift speed
6526a72c 2955 return fDriftSpeed[ fCluster.GetUniqueID() ];
2956}
2957
2958//________________________________________________________________________________________________________
2959Bool_t AliITSAlignMille2::FixedOrphans() const
2960{
2961 // are there fixed modules with no parent (normally in such a case
2962 // the constraints on the orphans should not be applied
2963 if (!IsConfigured()) {
2964 AliInfo("Still not configured");
2965 return kFALSE;
2966 }
2967 for (int i=0;i<fNModules;i++) {
2968 AliITSAlignMille2Module* md = GetMilleModule(i);
2969 if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
2970 }
2971 return kFALSE;
2972}
2973
2974//________________________________________________________________________________________________________
2975void AliITSAlignMille2::ConvertParamsToGlobal()
2976{
45993b23 2977 // convert params in local frame to global one
6526a72c 2978 double pars[AliITSAlignMille2Module::kMaxParGeom];
2979 for (int imd=fNModules;imd--;) {
2980 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2981 if (mod->GeomParamsGlobal()) continue;
2982 mod->GetGeomParamsGlo(pars);
2983 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
2984 mod->SetGeomParamsGlobal(kTRUE);
2985 }
2986}
2987
2988//________________________________________________________________________________________________________
2989void AliITSAlignMille2::ConvertParamsToLocal()
2990{
45993b23 2991 // convert params in global frame to local one
6526a72c 2992 double pars[AliITSAlignMille2Module::kMaxParGeom];
2993 for (int imd=fNModules;imd--;) {
2994 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2995 if (!mod->GeomParamsGlobal()) continue;
2996 mod->GetGeomParamsLoc(pars);
2997 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
2998 mod->SetGeomParamsGlobal(kFALSE);
2999 }
3000}
3001