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