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