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