]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSAlignMille2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Production file not written anymore
[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
ef24eb3b 25// fixing or constraning detection elements for best results.
6be22b3f 26//
27// author M. Lunardon (thanks to J. Castillo), ruben.shahoyan@cern.ch
28//-----------------------------------------------------------------------------
29
30#include <TFile.h>
ef24eb3b 31#include <TGrid.h>
6be22b3f 32#include <TClonesArray.h>
33#include <TMath.h>
34#include <TVirtualFitter.h>
35#include <TGeoManager.h>
36#include <TSystem.h>
37#include <TRandom.h>
38#include <TCollection.h>
39#include <TGeoPhysicalNode.h>
ef24eb3b 40#include <TMap.h>
41#include <TObjString.h>
42#include <TString.h>
6be22b3f 43#include "AliITSAlignMille2.h"
44#include "AliITSgeomTGeo.h"
45#include "AliGeomManager.h"
46#include "AliMillePede2.h"
47#include "AliTrackPointArray.h"
48#include "AliAlignObjParams.h"
49#include "AliLog.h"
50#include "AliTrackFitterRieman.h"
51#include "AliITSAlignMille2Constraint.h"
52#include "AliITSAlignMille2ConstrArray.h"
53#include "AliITSresponseSDD.h"
54#include "AliITSTPArrayFit.h"
55#include "AliCDBManager.h"
56#include "AliCDBStorage.h"
57#include "AliCDBEntry.h"
ef24eb3b 58#include "AliITSsegmentationSDD.h"
59#include "AliITSDriftSpeedArraySDD.h"
8102b2c9 60#include "AliITSCorrectSDDPoints.h"
ef24eb3b 61#include "AliESDVertex.h"
6be22b3f 62
63ClassImp(AliITSAlignMille2)
64
65const Char_t* AliITSAlignMille2::fgkRecKeys[] = {
66 "OCDB_PATH",
ef24eb3b 67 "OCDB_SPECIFIC",
6be22b3f 68 "GEOMETRY_FILE",
69 "SUPERMODULE_FILE",
70 "CONSTRAINTS_REFERENCE_FILE",
71 "PREALIGNMENT_FILE",
72 "PRECALIBSDD_FILE",
ef24eb3b 73 "PREVDRIFTSDD_FILE",
8102b2c9 74 "PRECORRMAPSDD_FILE",
75 "INITCORRMAPSDD_FILE",
6be22b3f 76 "INITCALBSDD_FILE",
ef24eb3b 77 "INITVDRIFTSDD_FILE",
6be22b3f 78 "INITDELTA_FILE",
8102b2c9 79 "INITGEOM_FILE",
6be22b3f 80 "SET_GLOBAL_DELTAS",
81 "CONSTRAINT_LOCAL",
82 "MODULE_VOLUID",
83 "MODULE_INDEX",
84 "SET_PSEUDO_PARENTS",
85 "SET_TRACK_FIT_METHOD",
86 "SET_MINPNT_TRA",
87 "SET_NSTDDEV",
88 "SET_RESCUT_INIT",
89 "SET_RESCUT_OTHER",
90 "SET_LOCALSIGMAFACTOR",
91 "SET_STARTFAC",
ef24eb3b 92 "SET_FINALFAC",
6be22b3f 93 "SET_B_FIELD",
94 "SET_SPARSE_MATRIX",
95 "REQUIRE_POINT",
96 "CONSTRAINT_ORPHANS",
97 "CONSTRAINT_SUBUNITS",
98 "APPLY_CONSTRAINT",
99 "SET_EXTRA_CLUSTERS_MODE",
100 "SET_USE_TPAFITTER",
101 "SET_USE_LOCAL_YERROR",
ef24eb3b 102 "SET_MIN_POINTS_PER_MODULE",
103 "SET_USE_SDDVDCORRMULT",
104 "SET_WEIGHT_PT",
105 "SET_USE_DIAMOND",
8102b2c9 106 "CORRECT_DIAMOND",
107 "SET_USE_VERTEX",
ef24eb3b 108 "SET_SAME_SDDT0"
6be22b3f 109};
110
111const Char_t AliITSAlignMille2::fgkXYZ[] = "XYZ";
112
113//========================================================================================================
114
115AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;
116Int_t AliITSAlignMille2::fgInstanceID = 0;
117
118//________________________________________________________________________________________________________
119AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename,TList *userInfo )
120: TObject(),
121 fMillepede(0),
122 fStartFac(16.),
ef24eb3b 123 fFinalFac(1.),
6be22b3f 124 fResCutInitial(100.),
125 fResCut(100.),
126 fNGlobal(0),
127 fNLocal(4),
128 fNStdDev(3),
129 fIsMilleInit(kFALSE),
130 fAllowPseudoParents(kFALSE),
131 //
132 fTPAFitter(0),
133 fCurrentModule(0),
134 fTrack(0),
135 fTrackBuff(0),
136 fCluster(),
137 fCurrentSensID(-1),
138 fClusLoc(12*3),
139 fClusGlo(12*3),
140 fClusSigLoc(12*3),
141 fGlobalDerivatives(0),
142 fMeasLoc(0),
143 fMeasGlo(0),
144 fSigmaLoc(0),
145 fConstrPT(-1),
146 fConstrPTErr(-1),
147 fConstrCharge(0),
8102b2c9 148 fRunID(0),
6be22b3f 149 //
150 fMinNPtsPerTrack(3),
ef24eb3b 151 fIniTrackParamsMeth(1),
6be22b3f 152 fTotBadLocEqPoints(0),
153 fRieman(0),
154 //
155 fConstraints(0),
156 fCacheMatrixOrig(kMaxITSSensID+1),
157 fCacheMatrixCurr(kMaxITSSensID+1),
158 //
159 fUseGlobalDelta(kFALSE),
6be22b3f 160 fTempExcludedModule(-1),
8102b2c9 161 fUserProvided(0),
6be22b3f 162 //
ef24eb3b 163 fIniUserInfo(userInfo),
8102b2c9 164 fIniGeomPath(""),
ef24eb3b 165 fIniDeltaPath(""),
166 fIniSDDRespPath(""),
6be22b3f 167 fPreCalSDDRespPath(""),
ef24eb3b 168 fIniSDDVDriftPath(""),
169 fPreSDDVDriftPath(""),
8102b2c9 170 fIniSDDCorrMapPath(""),
171 fPreSDDCorrMapPath(""),
172 fConvertPreDeltas(kFALSE),
ef24eb3b 173 fGeometryPath(""),
6be22b3f 174 fPreDeltaPath(""),
175 fConstrRefPath(""),
ef24eb3b 176 fDiamondPath(""),
6be22b3f 177 fGeoManager(0),
178 fIsConfigured(kFALSE),
179 fPreAlignQF(0),
180//
ef24eb3b 181 fIniRespSDD(0),
182 fPreRespSDD(0),
183 fIniVDriftSDD(0),
184 fPreVDriftSDD(0),
8102b2c9 185 fIniCorrMapSDD(0),
186 fPreCorrMapSDD(0),
ef24eb3b 187 fSegmentationSDD(0),
6be22b3f 188 fPrealignment(0),
189 fConstrRef(0),
190 fMilleModule(2),
191 fSuperModule(2),
192 fNModules(0),
193 fNSuperModules(0),
194 fUsePreAlignment(kFALSE),
195 fUseLocalYErr(kFALSE),
196 fBOn(kFALSE),
197 fBField(0.0),
1d06ac63 198 fDataType(kCosmics),
6be22b3f 199 fMinPntPerSens(0),
200 fBug(0),
201 fMilleVersion(2),
ef24eb3b 202 fExtraClustersMode(0),
203 fTrackWeight(1),
204 fWeightPt(0.),
205 fIsSDDVDriftMult(kFALSE),
206 fDiamond(),
207 fDiamondI(),
208 fUseDiamond(kFALSE),
8102b2c9 209 fUseVertex(kFALSE),
210 fVertexSet(kFALSE),
ef24eb3b 211 fDiamondPointID(-1),
8102b2c9 212 fDiamondModID(-1),
213 fCheckDiamondPoint(kDiamondCheckIfPrompt),
8693b6b7 214 fFixCurvIfConstraned(kTRUE),
215 fCurvFitWasConstrained(kFALSE),
8102b2c9 216 fConvAlgMatOld(100)
6be22b3f 217{
218 /// main constructor that takes input from configuration file
219 for (int i=3;i--;) fSigmaFactor[i] = 1.0;
220 //
221 // new RS
1d06ac63 222 for (int i=0;i<3;i++) {
8102b2c9 223 fCorrDiamond[i] = 0;
6be22b3f 224 }
1d06ac63 225 for (int itp=0;itp<kNDataType;itp++) {
226 fRequirePoints[itp] = kFALSE;
227 for (Int_t i=0; i<6; i++) {
228 fNReqLayUp[itp][i]=0;
229 fNReqLayDown[itp][i]=0;
230 fNReqLay[itp][i]=0;
231 }
232 for (Int_t i=0; i<3; i++) {
233 fNReqDetUp[itp][i]=0;
234 fNReqDetDown[itp][i]=0;
235 fNReqDet[itp][i]=0;
236 }
6be22b3f 237 }
238 //
ef24eb3b 239 // if (ProcessUserInfo(userInfo)) exit(1);
240 //
241 fDiamond.SetVolumeID(kVtxSensVID);
242 fDiamondI.SetVolumeID(kVtxSensVID);
243 float xyzd[3] = {0,0,0};
244 float covd[6] = {1,0,0,1,0,1e4};
245 fDiamond.SetXYZ(xyzd,covd); // dummy diamond
246 covd[5] = 1e-4;
247 fDiamondI.SetXYZ(xyzd,covd);
6be22b3f 248 //
249 Int_t lc=LoadConfig(configFilename);
250 if (lc) {
251 AliError(Form("Error %d loading configuration from %s",lc,configFilename));
252 exit(1);
253 }
254 //
255 fMillepede = new AliMillePede2();
256 fgInstance = this;
257 fgInstanceID++;
8102b2c9 258 ResetCovIScale();
6be22b3f 259 //
260}
261
262//________________________________________________________________________________________________________
263AliITSAlignMille2::~AliITSAlignMille2()
264{
265 /// Destructor
266 delete fMillepede;
267 delete[] fGlobalDerivatives;
268 delete fRieman;
269 delete fPrealignment;
270 delete fConstrRef;
ef24eb3b 271 delete fPreRespSDD;
272 delete fIniRespSDD;
273 delete fSegmentationSDD;
274 delete fIniVDriftSDD;
275 delete fPreVDriftSDD;
8102b2c9 276 delete fIniCorrMapSDD;
277 delete fPreCorrMapSDD;
6be22b3f 278 delete fTPAFitter;
279 fCacheMatrixOrig.Delete();
280 fCacheMatrixCurr.Delete();
281 fTrackBuff.Delete();
282 fConstraints.Delete();
283 fMilleModule.Delete();
284 fSuperModule.Delete();
285 if (--fgInstanceID==0) fgInstance = 0;
286}
287
288///////////////////////////////////////////////////////////////////////
289
290//________________________________________________________________________________________________________
291TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
292{
293 // read new record from config file
294 TString record;
295 static TObjArray* recElems = 0;
296 if (recElems) {delete recElems; recElems = 0;}
ef24eb3b 297 recOpt = "";
6be22b3f 298 //
299 TString keyws = recTitle;
300 if (!keyws.IsNull()) {
301 keyws.ToUpper();
302 // keyws += " ";
303 }
304 while (record.Gets(stream)) {
305 int cmt=record.Index("#");
306 if (cmt>=0) record.Remove(cmt); // skip comment
307 record.ReplaceAll("\t"," ");
308 record.ReplaceAll("\r"," ");
309 record.Remove(TString::kBoth,' ');
310 if (record.IsNull()) continue; // nothing to decode
311 if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
312 //
313 recElems = record.Tokenize(" ");
314 recTitle = recElems->At(0)->GetName();
315 recTitle.ToUpper();
316 recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
317 break;
318 }
319 if (rew || !recElems) rewind(stream);
320 return recElems;
321}
322
323//________________________________________________________________________________________________________
324Int_t AliITSAlignMille2::CheckConfigRecords(FILE* stream)
325{
b80c197e 326 // check the correctness of the record
6be22b3f 327 TString record,recTitle;
328 int lineCnt = 0;
329 rewind(stream);
330 while (record.Gets(stream)) {
331 int cmt=record.Index("#");
332 lineCnt++;
333 if (cmt>=0) record.Remove(cmt); // skip comment
334 record.ReplaceAll("\t"," ");
335 record.ReplaceAll("\r"," ");
336 record.Remove(TString::kBoth,' ');
337 if (record.IsNull()) continue; // nothing to decode
338 // extract keyword
339 int spc = record.Index(" ");
340 if (spc>0) recTitle = record(0,spc);
341 else recTitle = record;
342 recTitle.ToUpper();
343 Bool_t strOK = kFALSE;
344 for (int ik=kNKeyWords;ik--;) if (recTitle == fgkRecKeys[ik]) {strOK = kTRUE; break;}
345 if (strOK) continue;
346 //
347 AliError(Form("Unknown keyword %s at line %d",
348 recTitle.Data(),lineCnt));
349 return -1;
350 //
351 }
352 //
353 rewind(stream);
354 return 0;
355}
356
357
358//________________________________________________________________________________________________________
359Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
360{
361 // return 0 if success
362 // 1 if error in module index or voluid
363 //
ef24eb3b 364 AliInfo(Form("Loading MillePede2 configuration from %s",cfile));
365 AliCDBManager::Instance()->SetCacheFlag(kFALSE);
6be22b3f 366 FILE *pfc=fopen(cfile,"r");
367 if (!pfc) return -1;
368 //
369 TString record,recTitle,recOpt,recExt;
370 Int_t nrecElems,irec;
371 TObjArray *recArr=0;
372 //
373 fNModules = 0;
374 Bool_t stopped = kFALSE;
375 //
376 if (CheckConfigRecords(pfc)<0) return -1;
377 //
378 while(1) {
379 //
380 // ============= 1: we read some important records in predefined order ================
381 //
ef24eb3b 382 recTitle = fgkRecKeys[kOCDBDefaultPath];
383 if ( GetConfigRecord(pfc,recTitle,recOpt,1) && !recOpt.IsNull() ) {
384 AliInfo(Form("Configuration sets OCDB default storage to %s",recOpt.Data()));
385 AliCDBManager::Instance()->SetDefaultStorage( gSystem->ExpandPathName(recOpt.Data()) );
386 TObjString* objStr = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue("default");
387 if (!objStr) {stopped = kTRUE; break;}
388 objStr->SetUniqueID(1); // mark as user set
6be22b3f 389 }
390 //
ef24eb3b 391 if (fIniUserInfo && ProcessUserInfo(fIniUserInfo)) { AliError("Failed to process intial User Info"); stopped = kTRUE; break;}
6be22b3f 392 //
393 recTitle = fgkRecKeys[kGeomFile];
ef24eb3b 394 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fGeometryPath = gSystem->ExpandPathName(recOpt.Data());
8102b2c9 395 if ( LoadGeometry(fGeometryPath) ) { AliError("Failed to find/load target ideal Geometry"); stopped = kTRUE; break;}
6be22b3f 396 //
ef24eb3b 397 // Do we use new TrackPointArray fitter ?
398 recTitle = fgkRecKeys[kTPAFitter];
399 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) fTPAFitter = new AliITSTPArrayFit(kNLocal);
6be22b3f 400 //
401 recTitle = fgkRecKeys[kSuperModileFile];
402 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
403 recOpt.IsNull() ||
ef24eb3b 404 gSystem->ExpandPathName(recOpt) ||
6be22b3f 405 gSystem->AccessPathName(recOpt.Data()) ||
406 LoadSuperModuleFile(recOpt.Data()))
407 { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
408 //
6be22b3f 409 recTitle = fgkRecKeys[kConstrRefFile]; // LOCAL_CONSTRAINTS are defined wrt these deltas
410 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
411 if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
412 else {
413 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
414 if ( SetConstraintWrtRef(recOpt.Data()) )
415 { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
416 }
417 }
418 //
8102b2c9 419 recTitle = fgkRecKeys[kInitGeomFile];
420 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull() ) {
421 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
422 fIniGeomPath = recOpt;
423 gSystem->ExpandPathName(fIniGeomPath);
424 fUserProvided |= kSameInitGeomBit;
425 AliInfo(Form("Configuration sets Production Geometry to %s",fIniGeomPath.Data()));
426 }
427 if (fIniGeomPath.IsNull()) fIniGeomPath = fGeometryPath;
6be22b3f 428 //
abd7ef79 429 recTitle = fgkRecKeys[kInitDeltaFile];
ef24eb3b 430 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull() ) {
6be22b3f 431 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
ef24eb3b 432 fIniDeltaPath = recOpt;
433 gSystem->ExpandPathName(fIniDeltaPath);
8102b2c9 434 fUserProvided |= kSameInitDeltasBit;
ef24eb3b 435 AliInfo(Form("Configuration sets Production Deltas to %s",fIniDeltaPath.Data()));
6be22b3f 436 }
6be22b3f 437 //
abd7ef79 438 recTitle = fgkRecKeys[kPreDeltaFile];
ef24eb3b 439 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
440 if (!recOpt.IsNull()) {
441 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
442 fPreDeltaPath = recOpt;
443 gSystem->ExpandPathName(fPreDeltaPath);
444 }
445 else if (!fIniDeltaPath.IsNull()) {
8102b2c9 446 AliInfo("PreAlignment Deltas keyword is present but empty, will set to Init Deltas of the first tree");
ef24eb3b 447 fPreDeltaPath = fIniDeltaPath;
8102b2c9 448 if (fIniGeomPath != fGeometryPath) fConvertPreDeltas = kTRUE; // production and target geometries are different, request conversion
ef24eb3b 449 }
abd7ef79 450 AliInfo(Form("Configuration sets PreAlignment Deltas to %s",fPreDeltaPath.Data()));
6be22b3f 451 }
8102b2c9 452 //
453 // if initial deltas were provided, load them, apply to geometry and store are "original" matrices
454 if (CacheMatricesOrig()) {stopped = kTRUE; break;}
455 //
456 // then load prealignment deltas
457 if (!fPreDeltaPath.IsNull()) {
458 if (fConvertPreDeltas) ConvertDeltas(); // Prealignment deltas are the same as production ones, but need conversion to new geometry
459 else if (LoadDeltas(fPreDeltaPath,fPrealignment)) {stopped = kTRUE; break;} // read deltas from the file
460 }
abd7ef79 461 if (fPrealignment && ApplyToGeometry()) {stopped = kTRUE; break;}
6be22b3f 462 //
ef24eb3b 463 recTitle = fgkRecKeys[ kInitCalSDDFile ];
464 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
465 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
466 fIniSDDRespPath = recOpt;
467 gSystem->ExpandPathName(fIniSDDRespPath);
8102b2c9 468 fUserProvided |= kSameInitSDDRespBit;
ef24eb3b 469 AliInfo(Form("Configuration sets Production SDD Response to %s",fIniSDDRespPath.Data()));
470 }
471 if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {stopped = kTRUE; break;}
6be22b3f 472 //
8102b2c9 473 //
474 recTitle = fgkRecKeys[ kInitCorrMapSDDFile ];
475 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
476 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
477 fIniSDDCorrMapPath = recOpt;
478 gSystem->ExpandPathName(fIniSDDCorrMapPath);
479 fUserProvided |= kSameInitSDDCorrMapBit;
480 AliInfo(Form("Configuration sets Production SDD Correction Map to %s",fIniSDDCorrMapPath.Data()));
481 }
482 if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {stopped = kTRUE; break;}
483 //
6be22b3f 484 recTitle = fgkRecKeys[kPreCalSDDFile];
ef24eb3b 485 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
486 if (!recOpt.IsNull()) {
487 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
488 fPreCalSDDRespPath = recOpt;
489 gSystem->ExpandPathName(fPreCalSDDRespPath);
490 }
491 else if (!fIniSDDRespPath.IsNull()) {
492 AliInfo("PreCalibration SDD response keyword is present but empty, will set to Init SDD repsonse");
493 fPreCalSDDRespPath = fIniSDDRespPath;
494 }
6be22b3f 495 AliInfo(Form("Configuration sets PreCalibration SDD Response to %s",fPreCalSDDRespPath.Data()));
496 }
6be22b3f 497 //
ef24eb3b 498 if (LoadSDDResponse(fPreCalSDDRespPath, fPreRespSDD) ) {stopped = kTRUE; break;}
499 //
8102b2c9 500 recTitle = fgkRecKeys[kPreCorrMapSDDFile];
501 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) ) {
502 if (!recOpt.IsNull()) {
503 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
504 fPreSDDCorrMapPath = recOpt;
505 gSystem->ExpandPathName(fPreSDDCorrMapPath);
506 }
507 else if (!fIniSDDCorrMapPath.IsNull()) {
508 AliInfo("PreCalibration SDD Correction Map keyword is present but empty, will set to Init SDD Correction Map");
509 fPreSDDCorrMapPath = fIniSDDCorrMapPath;
510 }
511 AliInfo(Form("Configuration sets PreCalibration SDD Correction Map to %s",fPreSDDCorrMapPath.Data()));
512 }
ef24eb3b 513 //
8102b2c9 514 if (LoadSDDCorrMap(fPreSDDCorrMapPath, fPreCorrMapSDD) ) {stopped = kTRUE; break;}
515 // //
ef24eb3b 516 recTitle = fgkRecKeys[ kInitVDriftSDDFile ];
6be22b3f 517 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
518 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
ef24eb3b 519 fIniSDDVDriftPath = recOpt;
520 gSystem->ExpandPathName(fIniSDDVDriftPath);
8102b2c9 521 fUserProvided |= kSameInitSDDVDriftBit;
ef24eb3b 522 AliInfo(Form("Configuration sets Production SDD VDrift to %s",fIniSDDVDriftPath.Data()));
6be22b3f 523 }
ef24eb3b 524 if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {stopped = kTRUE; break;}
6be22b3f 525 //
ef24eb3b 526 recTitle = fgkRecKeys[ kPreVDriftSDDFile ];
527 if ( (recArr = GetConfigRecord(pfc,recTitle,recOpt,1)) && !recOpt.IsNull()) {
528 for (int i=2;i<=recArr->GetLast();i++) {recOpt += " "; recOpt += recArr->At(i)->GetName();} // in case of OCDB string
529 fPreSDDVDriftPath = recOpt;
530 gSystem->ExpandPathName(fPreSDDVDriftPath);
531 AliInfo(Form("Configuration sets PreCalibration SDD VDrift to %s",fPreSDDVDriftPath.Data()));
532 if (LoadSDDVDrift(fPreSDDVDriftPath, fPreVDriftSDD) ) {stopped = kTRUE; break;}
533 }
6be22b3f 534 //
535 recTitle = fgkRecKeys[ kGlobalDeltas ];
536 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
537 //
ef24eb3b 538 recTitle = fgkRecKeys[ kUseDiamond ];
539 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
540 if (!GetUseGlobalDelta()) {
541 AliError("Diamond constraint is supported only for Global Frame mode");
542 stopped = kTRUE;
543 break;
544 }
545 fUseDiamond = kTRUE;
546 if (!recOpt.IsNull()) {
547 fDiamondPath = recOpt;
548 gSystem->ExpandPathName(fDiamondPath);
8102b2c9 549 fUserProvided |= kSameDiamondBit;
ef24eb3b 550 AliInfo(Form("Configuration sets Diamond constraint to %s",fDiamondPath.Data()));
551 }
552 }
8102b2c9 553 //
554 recTitle = fgkRecKeys[ kUseVertex ];
555 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
556 if (!GetUseGlobalDelta()) {
557 AliError("Vertex constraint is supported only for Global Frame mode");
558 stopped = kTRUE;
559 break;
560 }
561 fUseVertex = kTRUE;
562 if (fUseDiamond) {
563 AliError("Cannot use Vertex and Diamond constraints at the same time");
564 stopped = kTRUE;
565 break;
566 }
567 AliInfo("Will use Vertex constraint when available");
568 }
6be22b3f 569 // =========== 2: see if there are local gaussian constraints defined =====================
570 // Note that they should be loaded before the modules declaration
571 //
572 recTitle = fgkRecKeys[ kConstrLocal ];
573 while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
574 nrecElems = recArr->GetLast()+1;
575 if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
576 if (GetConstraint(recOpt.Data())) {
577 AliError(Form("Existing constraint %s repeated",recOpt.Data()));
578 stopped = kTRUE; break;
579 }
580 recExt = recArr->At(2)->GetName();
581 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
582 double val = recExt.Atof();
583 recExt = recArr->At(3)->GetName();
584 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
585 double err = recExt.Atof();
586 int nwgh = nrecElems - 4;
587 double *wgh = new double[nwgh];
588 for (nwgh=0,irec=4;irec<nrecElems;irec++) {
589 recExt = recArr->At(irec)->GetName();
590 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
591 wgh[nwgh++] = recExt.Atof();
592 }
593 if (stopped) {delete[] wgh; break;}
594 //
595 ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
596 delete[] wgh;
597 //
598 } // end while for loop over local constraints
599 if (stopped) break;
600 //
601 // =========== 3: now read modules to align ===================================
602 //
603 rewind(pfc);
8fd71c0a 604 // create fixed modules
605 for (int j=0; j<fNSuperModules; j++) {
606 AliITSAlignMille2Module* proto = GetSuperModule(j);
607 if (!proto->IsAlignable()) continue;
608 AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(*proto);
609 // the matrix might be updated in case some prealignment was applied, check
610 TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
611 if (mup) *(mod->GetMatrix()) = *mup;
612 fMilleModule.AddAtAndExpand(mod,fNModules);
613 mod->SetGeomParamsGlobal(fUseGlobalDelta);
614 mod->SetUniqueID(fNModules++);
615 mod->SetNotInConf(kTRUE);
616 }
ef24eb3b 617 CreateVertexModule();
8fd71c0a 618 //
6be22b3f 619 while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
620 if (!(recTitle==fgkRecKeys[ kModVolID ] || recTitle==fgkRecKeys[ kModIndex ])) continue;
621 // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ] extra params]
622 // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
623 // sig* is the scaling parameters for the errors of the clusters of this module
624 // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
625 //
626 nrecElems = recArr->GetLast()+1;
627 if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
628 int idx = recOpt.Atoi();
629 UShort_t voluid = (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
630 AliITSAlignMille2Module* mod = 0;
631 //
632 if (voluid>=kMinITSSupeModuleID) { // custom supermodule
8fd71c0a 633 mod = GetMilleModuleByVID(voluid);
634 if (!mod) { // need to create
635 for (int j=0; j<fNSuperModules; j++) {
636 if (voluid==GetSuperModule(j)->GetVolumeID()) {
637 mod = new AliITSAlignMille2Module(*GetSuperModule(j));
638 // the matrix might be updated in case some prealignment was applied, check
639 TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
640 if (mup) *(mod->GetMatrix()) = *mup;
641 fMilleModule.AddAtAndExpand(mod,fNModules);
642 mod->SetGeomParamsGlobal(fUseGlobalDelta);
643 mod->SetUniqueID(fNModules++);
644 break;
645 }
646 }
6be22b3f 647 }
8fd71c0a 648 mod->SetNotInConf(kFALSE);
6be22b3f 649 }
650 else if (idx<=kMaxITSSensVID) {
651 mod = new AliITSAlignMille2Module(voluid);
652 fMilleModule.AddAtAndExpand(mod,fNModules);
8fd71c0a 653 mod->SetGeomParamsGlobal(fUseGlobalDelta);
654 mod->SetUniqueID(fNModules++);
6be22b3f 655 }
656 if (!mod) {stopped = kTRUE; break;} // bad volid
657 //
658 // geometry variation settings
659 for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
660 irec = i+2;
661 if (irec >= nrecElems) break;
662 recExt = recArr->At(irec)->GetName();
663 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
664 mod->SetFreeDOF(i, recExt.Atof() );
665 }
666 if (stopped) break;
667 //
668 // scaling factors for cluster errors
669 // first set default ones
670 for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);
671 for (int i=0;i<3;i++) {
672 irec = i+8;
673 if (irec >= nrecElems) break;
674 recExt = recArr->At(irec)->GetName();
675 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
676 mod->SetSigmaFactor(i, recExt.Atof() );
677 }
678 if (stopped) break;
679 //
6be22b3f 680 // now comes special detectors treatment
681 if (mod->IsSDD()) {
682 double vl = 0;
683 if (nrecElems>11) {
684 recExt = recArr->At(11)->GetName();
685 if (recExt.IsFloat()) vl = recExt.Atof();
686 else {stopped = kTRUE; break;}
687 irec = 11;
688 }
689 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
690 //
ef24eb3b 691 Bool_t cstLR = kFALSE;
692 for (int lr=0;lr<2;lr++) { // left right side vdrift corrections
693 vl = 0;
694 if (nrecElems>12+lr) {
695 recExt = recArr->At(12+lr)->GetName();
696 if (recExt.IsFloat()) vl = recExt.Atof();
697 else {stopped = kTRUE; break;}
698 irec = 12+lr;
699 }
700 mod->SetFreeDOF(lr==0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR,vl);
701 if (lr==1 && vl>=10) cstLR = kTRUE; // the right side should be constrained to left one
6be22b3f 702 }
ef24eb3b 703 if (cstLR) mod->SetVDriftLRSame();
6be22b3f 704 }
705 //
6be22b3f 706 mod->EvaluateDOF();
6be22b3f 707 //
708 // now check if there are local constraints on this module
709 for (++irec;irec<nrecElems;irec++) {
710 recExt = recArr->At(irec)->GetName();
711 if (recExt.IsFloat()) {stopped=kTRUE;break;}
712 AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
713 if (!cstr) {
714 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
715 stopped=kTRUE;
716 break;
717 }
718 cstr->AddModule(mod);
719 }
720 if (stopped) break;
721 } // end while for loop over modules
722 if (stopped) break;
723 //
724 if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}
725 BuildHierarchy(); // preprocess loaded modules
726 //
727 // =========== 4: the rest may come in arbitrary order =======================================
728 rewind(pfc);
729 while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
730 //
731 nrecElems = recArr->GetLast()+1;
732 //
733 // some simple flags -----------------------------------------------------------------------
734 //
735 if (recTitle == fgkRecKeys[ kPseudoParents ]) SetAllowPseudoParents(kTRUE);
736 //
737 // some optional parameters ----------------------------------------------------------------
738 else if (recTitle == fgkRecKeys[ kTrackFitMethod ]) {
739 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
740 SetInitTrackParamsMeth(recOpt.Atoi());
741 }
742 //
743 else if (recTitle == fgkRecKeys[ kMinPntTrack ]) {
744 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
745 fMinNPtsPerTrack = recOpt.Atoi();
746 }
747 //
748 else if (recTitle == fgkRecKeys[ kNStDev ]) {
749 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
750 fNStdDev = (Int_t)recOpt.Atof();
751 }
752 //
753 else if (recTitle == fgkRecKeys[ kResCutInit ]) {
754 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
755 fResCutInitial = recOpt.Atof();
756 }
757 //
758 else if (recTitle == fgkRecKeys[ kResCutOther ]) {
759 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
760 fResCut = recOpt.Atof();
761 }
762 //
763 else if (recTitle == fgkRecKeys[ kLocalSigFactor ]) { //-------------------------
764 for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
765 fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
766 if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
767 }
768 if (stopped) break;
769 }
770 //
771 else if (recTitle == fgkRecKeys[ kStartFactor ]) { //-------------------------
772 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
773 fStartFac = recOpt.Atof();
774 }
775 //
ef24eb3b 776 else if (recTitle == fgkRecKeys[ kFinalFactor ]) { //-------------------------
777 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
778 fFinalFac = recOpt.Atof();
779 }
780 //
6be22b3f 781 // pepo2708909
782 else if (recTitle == fgkRecKeys[ kExtraClustersMode ]) { //-------------------------
783 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
784 fExtraClustersMode = recOpt.Atoi();
785 }
786 // endpepo270809
787 //
788 else if (recTitle == fgkRecKeys[ kBField ]) { //-------------------------
789 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
790 SetBField( recOpt.Atof() );
791 }
792 //
ef24eb3b 793 else if (recTitle == fgkRecKeys[ kSDDVDCorrMult ]) { //-------------------------
794 SetSDDVDCorrMult( recOpt.IsNull() || (recOpt.IsFloat() && (recOpt.Atof())>-0.5) );
795 }
796 //
797 else if (recTitle == fgkRecKeys[ kWeightPt ]) { //-------------------------
798 double wgh = 1;
799 if (!recOpt.IsNull()) {
800 if (!recOpt.IsFloat()) {stopped = kTRUE; break;}
801 else wgh = recOpt.Atof();
802 }
803 SetWeightPt(wgh);
804 }
805 //
6be22b3f 806 else if (recTitle == fgkRecKeys[ kSparseMatrix ]) { // matrix solver type
807 //
808 AliMillePede2::SetGlobalMatSparse(kTRUE);
809 if (recOpt.IsNull()) continue;
810 // solver type and settings
811 if (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
812 else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
813 else {stopped = kTRUE; break;}
814 //
815 if (nrecElems>=3) { // preconditioner type
816 recExt = recArr->At(2)->GetName();
817 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
818 AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
819 }
820 //
821 if (nrecElems>=4) { // tolerance
822 recExt = recArr->At(3)->GetName();
823 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
824 AliMillePede2::SetMinResTol( recExt.Atof() );
825 }
826 //
827 if (nrecElems>=5) { // maxIter
828 recExt = recArr->At(4)->GetName();
829 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
830 AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
831 }
832 }
833 //
834 else if (recTitle == fgkRecKeys[ kRequirePoint ]) { //-------------------------
835 // syntax: REQUIRE_POINT where ndet updw nreqpts
836 // where = LAYER or DETECTOR
837 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
838 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
839 // nreqpts = minimum number of points of that type
840 if (nrecElems>=5) {
841 recOpt.ToUpper();
842 int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
843 int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
844 int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
1d06ac63 845 //
846 int rtp = -1; // use for run type
847 if (nrecElems>5) {
848 TString tpstr = ((TObjString*)recArr->At(5))->GetString();
849 if ( tpstr.Contains("cosmics",TString::kIgnoreCase) ) rtp = kCosmics;
850 else if ( tpstr.Contains("collision",TString::kIgnoreCase) ) rtp = kCollision;
851 else {stopped = kTRUE; break;}
6be22b3f 852 }
1d06ac63 853 //
854 int tpmn= rtp<0 ? 0 : rtp;
855 int tpmx= rtp<0 ? kNDataType-1 : rtp;
856 for (int itp=tpmn;itp<=tpmx;itp++) {
857 fRequirePoints[itp]=kTRUE;
858 if (recOpt == "LAYER") {
859 if (lr<0 || lr>5) {stopped = kTRUE; break;}
860 if (hb>0) fNReqLayUp[itp][lr]=np;
861 else if (hb<0) fNReqLayDown[itp][lr]=np;
862 else fNReqLay[itp][lr]=np;
863 }
864 else if (recOpt == "DETECTOR") {
865 if (lr<0 || lr>2) {stopped = kTRUE; break;}
866 if (hb>0) fNReqDetUp[itp][lr]=np;
867 else if (hb<0) fNReqDetDown[itp][lr]=np;
868 else fNReqDet[itp][lr]=np;
869 }
870 else {stopped = kTRUE; break;}
6be22b3f 871 }
1d06ac63 872 if (stopped) break;
6be22b3f 873 }
874 else {stopped = kTRUE; break;}
875 }
876 //
877 // global constraints on the subunits/orphans
878 else if (recTitle == fgkRecKeys[ kConstrOrphans ]) { //------------------------
879 // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
880 if (nrecElems<4) {stopped = kTRUE; break;}
881 recExt = recArr->At(2)->GetName();
882 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
883 double val = recExt.Atof();
884 UInt_t pattern = 0;
885 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
886 recExt = recArr->At(irec)->GetName();
887 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
888 pattern |= 0x1 << recExt.Atoi();
889 }
890 if (stopped) break;
891 if (recOpt == "MEAN") ConstrainOrphansMean(val,pattern);
892 else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
893 else {stopped = kTRUE; break;}
894 }
895 //
896 else if (recTitle == fgkRecKeys[ kConstrSubunits ]) { //------------------------
ef24eb3b 897 // expect CONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
6be22b3f 898 if (nrecElems<5) {stopped = kTRUE; break;}
899 recExt = recArr->At(2)->GetName();
900 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
901 double val = recExt.Atof();
902 UInt_t pattern = 0;
903 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
904 recExt = recArr->At(irec)->GetName();
905 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
906 int parid = recExt.Atoi();
907 if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
908 else break; // list of params is over
909 }
910 if (stopped) break;
911 //
912 Bool_t meanC;
913 if (recOpt == "MEAN") meanC = kTRUE;
914 else if (recOpt == "MEDIAN") meanC = kFALSE;
915 else {stopped = kTRUE; break;}
916 //
917 int curID = -1;
918 int rangeStart = -1;
919 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
920 recExt = recArr->At(irec)->GetName();
921 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
922 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
923 else curID = recExt.Atoi();
924 //
925 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
926 // this was a range start or single
927 int start;
928 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
929 else start = curID; // create constraint either for single module (or 1st in the range)
930 for (int id=start;id<=curID;id++) {
931 int id0 = IsVIDDefined(id);
932 if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
933 if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
934 else ConstrainModuleSubUnitsMedian(id0,val,pattern);
935 }
936 }
937 if (rangeStart>=0) stopped = kTRUE; // unfinished range
938 if (stopped) break;
939 }
940 //
941 // association of modules with local constraints
942 else if (recTitle == fgkRecKeys[ kApplyConstr ]) { //------------------------
943 // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
944 if (nrecElems<3) {stopped = kTRUE; break;}
945 int nmID0=-1,nmID1=-1;
946 for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
947 recExt = recArr->At(irec)->GetName();
948 if (recExt.IsFloat()) break;
949 // check if such a constraint was declared
950 if (!GetConstraint(recExt.Data())) {
951 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
952 stopped=kTRUE;
953 break;
954 }
955 if (nmID0<0) nmID0 = irec;
956 nmID1 = irec;
957 }
958 if (stopped) break;
959 //
960 if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
961 //
962 // now read the list of modules to constrain
963 int curID = -1;
964 int rangeStart = -1;
965 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
966 recExt = recArr->At(irec)->GetName();
967 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
968 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
969 else curID = recExt.Atoi();
970 //
971 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
972 //
973 // this was a range start or single
974 int start;
975 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
976 else start = curID; // create constraint either for single module (or 1st in the range)
977 for (int id=start;id<=curID;id++) {
978 AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
979 if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
980 for (int nmid=nmID0;nmid<=nmID1;nmid++)
981 ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
982 }
983 }
984 if (rangeStart>=0) stopped = kTRUE; // unfinished range
985 if (stopped) break;
986 }
ef24eb3b 987 //
988 // request of the same T0 for group of SDD modules
989 else if (recTitle == fgkRecKeys[ kSameSDDT0 ]) { //------------------------
990 // expect SET_SAME_SDDT0 [SensID1 ... SensIDn - SensIDm]
991 if (nrecElems<3) {stopped = kTRUE; break;}
992 //
993 // now read the list of modules to constrain
994 int curID = -1;
995 int rangeStart = -1;
996 AliITSAlignMille2ConstrArray *cstrT0 = new AliITSAlignMille2ConstrArray("SDDT0",0,0,0,0);
997 int naddM = 0;
998 cstrT0->SetPattern(BIT(AliITSAlignMille2Module::kDOFT0));
999 for (irec=1;irec<nrecElems;irec++) { // read modules to apply this constraint
1000 recExt = recArr->At(irec)->GetName();
1001 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
1002 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
1003 else curID = recExt.Atoi();
1004 //
1005 if (curID<kSDDoffsID || curID>=kSDDoffsID+kNSDDmod) {stopped = kTRUE; break;}
1006 //
1007 // this was a range start or single
1008 int start;
1009 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
1010 else start = curID; // create constraint either for single module (or 1st in the range)
1011 for (int id=start;id<=curID;id++) {
1012 int vid = AliITSAlignMille2Module::GetVolumeIDFromIndex(id);
1013 if (vid<=1) {AliDebug(3,Form("Undefined module index %d requested in the SAME_SDDT0 constraint, skipping",id)); continue;}
1014 AliITSAlignMille2Module *md = GetMilleModuleByVID(vid);
1015 if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
1016 cstrT0->AddModule(md,kFALSE);
1017 naddM++;
1018 }
1019 }
1020 if (rangeStart>=0) stopped = kTRUE; // unfinished range
1021 if (stopped) break;
1022 if (naddM<2) delete cstrT0;
1023 else {
1024 cstrT0->SetConstraintID(GetNConstraints());
1025 fConstraints.Add(cstrT0);
1026 }
6be22b3f 1027 }
ef24eb3b 1028 //
6be22b3f 1029 // Do we use new local Y errors?
1030 else if (recTitle == fgkRecKeys[ kUseLocalYErr ]) {
1031 // expect SET_TPAFITTER
1032 fUseLocalYErr = kTRUE;
1033 }
1034 //
1035 else if (recTitle == fgkRecKeys[ kMinPointsSens ]) { //-------------------------
1036 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
1037 SetMinPointsPerSensor( recOpt.Atoi() );
1038 }
1039 //
ef24eb3b 1040 else if (recTitle == fgkRecKeys[ kOCDBSpecificPath ]) { //-------------------------
1041 if (recOpt.IsNull() || nrecElems<3 ) {stopped = kTRUE; break;}
1042 AliCDBManager::Instance()->SetSpecificStorage(recOpt.Data(), gSystem->ExpandPathName(recArr->At(2)->GetName()));
1043 AliInfo(Form("Configuration sets OCDB specific storage %s to %s",recOpt.Data(),recArr->At(2)->GetName()));
1044 TObjString *pths = (TObjString*)AliCDBManager::Instance()->GetStorageMap()->GetValue(recOpt.Data());
1045 if (!pths) { stopped = kTRUE; break; }
1046 pths->SetUniqueID(1); // mark as set by user
1047 }
1048 //
8102b2c9 1049 else if (recTitle == fgkRecKeys[ kCorrectDiamond ] && fUseDiamond) { //-------------------------
1050 if (nrecElems<4) {stopped = kTRUE; break;}
1051 for (int i=0;i<3;i++) fCorrDiamond[i] = ((TObjString*)recArr->At(i+1))->GetString().Atof();
1052 AliInfo(Form("Correction %+.4f %+.4f %+.4f will be applied to diamond",fCorrDiamond[0],fCorrDiamond[1],fCorrDiamond[2]));
1053 }
1054 //
6be22b3f 1055 else continue; // already processed record
1056 //
1057 } // end of while loop 4 over the various params
1058 //
1059 break;
1060 } // end of while(1) loop
1061 //
1062 fclose(pfc);
ef24eb3b 1063 if (!fDiamondPath.IsNull() && IsDiamondUsed() && LoadDiamond(fDiamondPath) ) stopped = kTRUE;
6be22b3f 1064 if (stopped) {
1065 AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
1066 return -1;
1067 }
1068 //
abd7ef79 1069 if (CacheMatricesCurr()) return -1;
6be22b3f 1070 SetUseLocalYErrors(fUseLocalYErr); // YErr used only with TPAFitter
ef24eb3b 1071 fSegmentationSDD = new AliITSsegmentationSDD();
1072 //
6be22b3f 1073 fIsConfigured = kTRUE;
1074 return 0;
1075}
1076
1077//________________________________________________________________________________________________________
1078void AliITSAlignMille2::BuildHierarchy()
1079{
1080 // build the hieararhy of the modules to align
1081 //
1082 if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
1083 AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
1084 "Since Deltas are local, switching to NoPseudoParents");
1085 SetAllowPseudoParents(kFALSE);
1086 }
1087 // set parent/child relationship for modules to align
1088 AliInfo("Setting parent/child relationships\n");
1089 //
1090 // 1) child -> parent reference
1091 for (int ipar=0;ipar<fNModules;ipar++) {
1092 AliITSAlignMille2Module* parent = GetMilleModule(ipar);
1093 if (parent->IsSensor()) continue; // sensor cannot be a parent
1094 //
1095 for (int icld=0;icld<fNModules;icld++) {
1096 if (icld==ipar) continue;
1097 AliITSAlignMille2Module* child = GetMilleModule(icld);
1098 if (!child->BelongsTo(parent)) continue;
1099 // child cannot have more sensors than the parent
1100 if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
1101 //
1102 AliITSAlignMille2Module* parOld = child->GetParent();
1103 // is this parent candidate closer than the old parent ?
1104 if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
1105 child->SetParent(parent);
1106 }
1107 //
1108 }
1109 //
1110 // add parent -> children reference
1111 for (int icld=0;icld<fNModules;icld++) {
1112 AliITSAlignMille2Module* child = GetMilleModule(icld);
1113 AliITSAlignMille2Module* parent = child->GetParent();
1114 if (parent) parent->AddChild(child);
1115 }
1116 //
1117 // reorder the modules in such a way that parents come first
1118 for (int icld=0;icld<fNModules;icld++) {
1119 AliITSAlignMille2Module* child = GetMilleModule(icld);
1120 AliITSAlignMille2Module* parent;
1121 while ( (parent=child->GetParent()) && (parent->GetUniqueID()>child->GetUniqueID()) ) {
1122 // swap
1123 fMilleModule[icld] = parent;
1124 fMilleModule[parent->GetUniqueID()] = child;
1125 child->SetUniqueID(parent->GetUniqueID());
1126 parent->SetUniqueID(icld);
1127 child = parent;
1128 }
1129 //
1130 }
1131 //
1132 // Go over the child->parent chain and mark modules with explicitly provided sensors.
1133 // If the sensors of the unit are explicitly declared, all undeclared sensors are
1134 // suppresed in this unit.
1135 for (int icld=fNModules;icld--;) {
1136 AliITSAlignMille2Module* child = GetMilleModule(icld);
1137 AliITSAlignMille2Module* parent = child->GetParent();
1138 if (!parent) continue;
1139 //
1140 // check if this parent was already processed
1141 if (!parent->AreSensorsProvided()) {
1142 parent->DelSensitiveVolumes();
1143 parent->SetSensorsProvided(kTRUE);
1144 }
1145 // reattach sensors to parent
1146 for (int isc=child->GetNSensitiveVolumes();isc--;) {
1147 UShort_t senVID = child->GetSensVolVolumeID(isc);
1148 if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
1149 }
1150 }
1151 //
1152}
1153
1154// pepo
1155//________________________________________________________________________________________________________
1156void AliITSAlignMille2::SetCurrentModule(Int_t id)
1157{
1158 // set the current supermodule
1159 // new meaning
1160 if (fMilleVersion>=2) {
1161 fCurrentModule = GetMilleModule(id);
1162 return;
1163 }
1164 // old meaning
1165 if (fMilleVersion<=1) {
1166 Int_t index=id;
1167 /// set as current the SuperModule that contains the 'index' sens.vol.
1168 if (index<0 || index>2197) {
1169 AliInfo("index does not correspond to a sensitive volume!");
1170 return;
1171 }
1172 UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
1173 Int_t k=IsContained(voluid);
1174 if (k>=0){
1175 fCurrentSensID = index;
1176 fCluster.SetVolumeID(voluid);
1177 fCluster.SetXYZ(0,0,0);
1178 InitModuleParams();
1179 }
1180 else
1181 AliInfo(Form("module %d not defined\n",index));
1182 }
1183}
1184// endpepo
1185//________________________________________________________________________________________________________
1d06ac63 1186void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts,Int_t runtype)
6be22b3f 1187{
1188 // set minimum number of points in specific detector or layer
1189 // where = LAYER or DETECTOR
1190 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
1191 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
1192 // nreqpts = minimum number of points of that type
1193 ndet--;
1d06ac63 1194 int tpmn= runtype<0 ? 0 : runtype;
1195 int tpmx= runtype<0 ? kNDataType-1 : runtype;
1196 //
1197 for (int itp=tpmn;itp<=tpmx;itp++) {
1198 fRequirePoints[itp]=kTRUE;
1199 if (strstr(where,"LAYER")) {
1200 if (ndet<0 || ndet>5) return;
1201 if (updw>0) fNReqLayUp[itp][ndet]=nreqpts;
1202 else if (updw<0) fNReqLayDown[itp][ndet]=nreqpts;
1203 else fNReqLay[itp][ndet]=nreqpts;
1204 }
1205 else if (strstr(where,"DETECTOR")) {
1206 if (ndet<0 || ndet>2) return;
1207 if (updw>0) fNReqDetUp[itp][ndet]=nreqpts;
1208 else if (updw<0) fNReqDetDown[itp][ndet]=nreqpts;
1209 else fNReqDet[itp][ndet]=nreqpts;
1210 }
6be22b3f 1211 }
1212}
1213
1214//________________________________________________________________________________________________________
1215Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname)
1216{
1217 /// index from symname
1218 if (!symname) return -1;
1219 for (Int_t i=0;i<=kMaxITSSensID; i++) {
1220 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
1221 }
1222 return -1;
1223}
1224
1225//________________________________________________________________________________________________________
1226Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid)
1227{
1228 /// index from volume ID
1229 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
1230 if (lay<1|| lay>6) return -1;
1231 Int_t idx=Int_t(voluid)-2048*lay;
1232 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
1233 for (Int_t ilay=1; ilay<lay; ilay++)
1234 idx += AliGeomManager::LayerSize(ilay);
1235 return idx;
1236}
1237
1238//________________________________________________________________________________________________________
1239UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname)
1240{
1241 /// volume ID from symname
1242 /// works for sensitive volumes only
1243 if (!symname) return 0;
1244
1245 for (UShort_t voluid=2000; voluid<13300; voluid++) {
1246 Int_t modId;
1247 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
1248 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
1249 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
1250 }
1251 }
1252
1253 return 0;
1254}
1255
1256//________________________________________________________________________________________________________
1257UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index)
1258{
1259 /// volume ID from index
1260 if (index<0) return 0;
1261 if (index<2198)
1262 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
1263 else {
1264 for (int i=0; i<fNSuperModules; i++) {
1265 if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
1266 }
1267 }
1268 return 0;
1269}
1270
1271//________________________________________________________________________________________________________
8102b2c9 1272Int_t AliITSAlignMille2::LoadGeometry(TString& path)
6be22b3f 1273{
8102b2c9 1274 // initialize ideal geometry
1275 AliInfo(Form("Loading ideal geometry %s",path.Data()));
1276 if (path.IsNull()) {
1277 AliError("Path to geometry is not provided");
6be22b3f 1278 return -1;
1279 }
1280 //
8102b2c9 1281 AliCDBEntry *entry = 0;
1282 TGeoManager *gm = 0;
1283 while(1) {
1284 if (path.BeginsWith("path: ")) { // must load from OCDB
1285 entry = GetCDBEntry(path.Data());
1286 if (!entry) break;
1287 gm = (TGeoManager*) entry->GetObject();
1288 entry->SetObject(NULL);
1289 entry->SetOwner(kTRUE);
1290 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
1291 // delete cdbId;
1292 // delete entry;
1293 break;
1294 }
1295 //
1296 if (gSystem->AccessPathName(path.Data())) break;
1297 TFile* precf = TFile::Open(path.Data());
1298 if (precf->FindKey("ALICE")) gm = (TGeoManager*)precf->Get("ALICE");
1299 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
1300 gm = (TGeoManager*) entry->GetObject();
1301 if (gm && gm->InheritsFrom(TGeoManager::Class())) entry->SetObject(NULL);
1302 else gm = 0;
1303 entry->SetObject(NULL);
1304 entry->SetOwner(kTRUE);
1305 delete entry;
1306 }
1307 precf->Close();
1308 delete precf;
1309 break;
1310 }
1311 //
1312 if (!gm) {AliError(Form("Failed to load geometry from %s",path.Data())); return -1;}
1313 AliGeomManager::SetGeometry(gm);
6be22b3f 1314 fGeoManager = AliGeomManager::GetGeometry();
1315 if (!fGeoManager) {
1316 AliInfo("Couldn't initialize geometry");
1317 return -1;
1318 }
1319 return 0;
1320}
1321
1322//________________________________________________________________________________________________________
1323Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname)
1324{
1325 // Load the global deltas from this file. The local gaussian constraints on some modules
1326 // will be defined with respect to the deltas from this reference file, converted to local
1327 // delta format. Note: conversion to local format requires reloading the geometry!
1328 //
1329 AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
1330 if (!fGeoManager) return -1;
1331 fConstrRefPath = reffname;
1332 if (fConstrRefPath == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
1333 fConstrRef = new TClonesArray("AliAlignObjParams",1);
1334 return 0;
1335 }
1336 if (LoadDeltas(fConstrRefPath,fConstrRef)) return -1;
1337 //
1338 // we need ideal geometry to convert global deltas to local ones
1339 if (fUsePreAlignment) {
1340 AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
1341 return -1;
1342 }
1343 //
1344 AliInfo("Converting global reference deltas to local ones");
1345 Int_t nprea = fConstrRef->GetEntriesFast();
1346 for (int ix=0; ix<nprea; ix++) {
1347 AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
1348 if (!preo->ApplyToGeometry()) return -1;
1349 }
1350 //
1351 // now convert the global reference deltas to local ones
1352 for (int i=fConstrRef->GetEntriesFast();i--;) {
1353 AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
1354 TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
1355 if (!mupd) { // this is not alignable entry, need to look in the supermodules
1356 for (int im=fNSuperModules;im--;) {
1357 AliITSAlignMille2Module* mod = GetSuperModule(im);
1358 if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
1359 mupd = mod->GetMatrix();
1360 break;
1361 }
1362 if (!mupd) {
1363 AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
1364 return -1;
1365 }
1366 }
1367 TGeoHMatrix preMat;
1368 preo->GetMatrix(preMat); // Delta_Glob
1369 TGeoHMatrix tmpMat = *mupd; // Delta_Glob * Delta_Glob_Par * M
1370 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
1371 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
1372 preo->SetMatrix(tmpMat); // local corrections
1373 }
1374 //
1375 // we need to reload the geometry spoiled by this reference deltas...
1376 delete fGeoManager;
8102b2c9 1377 AliInfo("Reloading target ideal geometry");
1378 return LoadGeometry(fGeometryPath);
6be22b3f 1379 //
1380}
1381
1382//________________________________________________________________________________________________________
1383void AliITSAlignMille2::Init()
1384{
1385 // perform global initialization
1386 //
1387 if (fIsMilleInit) {
1388 AliInfo("Millepede has been already initialized!");
1389 return;
1390 }
1391 // range constraints in such a way that the childs are constrained before their parents
1392 // orphan constraints come last
1393 for (int ic=0;ic<GetNConstraints();ic++) {
1394 for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
1395 AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
1396 AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
1397 if (cst0->GetModuleID()<cst1->GetModuleID()) {
1398 // swap
1399 fConstraints[ic] = cst1;
1400 fConstraints[ic1] = cst0;
1401 }
1402 }
1403 }
1404 //
1405 if (!GetUseGlobalDelta()) {
1406 AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
1407 for (int imd=fNModules;imd--;) {
1408 AliITSAlignMille2Module* mod = GetMilleModule(imd);
1409 int npar = mod->GetNParTot();
1410 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1411 for (int ipar=0;ipar<npar;ipar++) {
1412 if (!mod->IsFreeDOF(ipar)) continue;
1413 mod->SetParOffset(ipar,fNGlobal++);
1414 }
1415 }
1416 }
1417 else {
1418 // init millepede, decide which parameters are to be fitted explicitly
1419 for (int imd=fNModules;imd--;) {
1420 AliITSAlignMille2Module* mod = GetMilleModule(imd);
8102b2c9 1421 if (mod->IsNotInConf()) continue; // dummy module
6be22b3f 1422 int npar = mod->GetNParTot();
1423 // the parameter may have max 1 free instance, otherwise the equations are underdefined
1424 for (int ipar=0;ipar<npar;ipar++) {
1425 if (!mod->IsFreeDOF(ipar)) continue; // fixed
1426 //
1427 int nFreeInstances = 0;
1428 //
1429 AliITSAlignMille2Module* parent = mod;
1430 Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
1431 //
1432 Bool_t addToFit = kFALSE;
1433 // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
1434 // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
1435 // 2) the same applies to all of its parents
1436 // 3) it has at least 1 unconstrained direct child
1437 while(parent) {
8102b2c9 1438 if (parent->IsNotInConf()) {parent = parent->GetParent(); continue;}
6be22b3f 1439 if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
1440 nFreeInstances++;
1441 if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
1442 if (cstGauss) addToFit = kTRUE;
1443 parent = parent->GetParent();
1444 }
1445 if (nFreeInstances>1) {
1446 AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
1447 "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
1448 exit(1);
1449 }
1450 //
1451 // i) Are PseudoParents allowed?
1452 if (!PseudoParentsAllowed()) addToFit = kTRUE;
1453 // ii) check if this module has no child with such a free parameter. Since the order of this check
1454 // goes from child to parent, by this moment such a parameter must have been already added
1455 else if (!IsParModFamilyVaried(mod,ipar)) addToFit = kTRUE; // no varied children at all
1456 else if (!IsParFamilyFree(mod,ipar,1)) addToFit = kTRUE; // no unconstrained direct children
1457 // otherwise the value of this parameter can be extracted from simple contraint and the values of
1458 // the relevant parameters of its children the fit is done. Hence it is not included
1459 if (!addToFit) continue;
1460 //
1461 // shall add this parameter to explicit fit
1462 // printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
1463 mod->SetParOffset(ipar,fNGlobal++);
1464 }
1465 }
1466 }
1467 //
1468 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, kNLocal, fNStdDev));
1469 fGlobalDerivatives = new Double_t[fNGlobal];
1470 memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1471 //
1472 fMillepede->InitMille(fNGlobal,kNLocal,fNStdDev,fResCut,fResCutInitial);
1473 fMillepede->SetMinPntValid(fMinPntPerSens);
1474 fIsMilleInit = kTRUE;
1475 //
1476 ResetLocalEquation();
1477 AliInfo("Parameters initialized to zero");
1478 //
1479 /// Fix non free parameters
1480 for (Int_t i=0; i<fNModules; i++) {
1481 AliITSAlignMille2Module* mod = GetMilleModule(i);
1482 for (Int_t j=0; j<mod->GetNParTot(); j++) {
1483 if (mod->GetParOffset(j)<0) continue; // not varied
1484 FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1485 fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1486 }
1487 }
1488 //
8102b2c9 1489 ResetCovIScale();
6be22b3f 1490 // Set iterations
1491 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
ef24eb3b 1492 if (fFinalFac>1) fMillepede->SetChi2CutRef(fFinalFac);
8102b2c9 1493 fTrackBuff.Expand(24);
6be22b3f 1494 //
1495}
1496
1497//________________________________________________________________________________________________________
1498void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma)
1499{
1500 /// Constrain equation defined by par to value
1501 if (!fIsMilleInit) Init();
1502 fMillepede->SetGlobalConstraint(par, value, sigma);
1503 AliInfo("Adding constraint");
1504}
1505
1506//________________________________________________________________________________________________________
1507void AliITSAlignMille2::InitGlobalParameters(Double_t *par)
1508{
1509 /// Initialize global parameters with par array
1510 if (!fIsMilleInit) Init();
1511 fMillepede->SetGlobalParameters(par);
1512 AliInfo("Init Global Parameters");
1513}
1514
1515//________________________________________________________________________________________________________
1516void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value)
1517{
1518 /// Parameter iPar is encourage to vary in [-value;value].
1519 /// If value == 0, parameter is fixed
1520 if (!fIsMilleInit) {
1521 AliInfo("Millepede has not been initialized!");
1522 return;
1523 }
1524 fMillepede->SetParSigma(iPar, value);
1525 if (IsZero(value)) AliInfo(Form("Parameter %i Fixed", iPar));
1526}
1527
1528//________________________________________________________________________________________________________
1529void AliITSAlignMille2::ResetLocalEquation()
1530{
1531 /// Reset the derivative vectors
1532 for(int i=kNLocal;i--;) fLocalDerivatives[i] = 0.0;
1533 memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
1534}
1535
1536//________________________________________________________________________________________________________
1537Int_t AliITSAlignMille2::ApplyToGeometry()
1538{
1539 // apply prealignment to ideal geometry
1540 Int_t nprea = fPrealignment->GetEntriesFast();
1541 AliInfo(Form("Array of prealignment deltas: %d entries",nprea));
1542 //
1543 for (int ix=0; ix<nprea; ix++) {
1544 AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1545 Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1546 if (index>=0) {
1547 if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1548 fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1549 }
1550 if (!preo->ApplyToGeometry()) {
1551 AliError(Form("Failed on ApplyToGeometry at %s",preo->GetSymName()));
1552 return -1;
1553 }
1554 }
1555 //
1556 fUsePreAlignment = kTRUE;
1557 return 0;
1558}
1559
1560//________________________________________________________________________________________________________
1561Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1562{
1563 // quality factors from prealignment
1564 if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1565 return fPreAlignQF[index]-1;
1566}
1567
1568//________________________________________________________________________________________________________
1569AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp)
1570{
1571 /// create a new AliTrackPointArray keeping only defined modules
1572 /// move points according to a given prealignment, if any
8102b2c9 1573 /// sort alitrackpoints w.r.t. global Y direction, if cosmics
6be22b3f 1574 const Double_t kRad2L[6] = {5*5,10*10,18*18,30*30,40*40,60*60};
1575 const Float_t kSensSigY2[6] = {200e-4*200e-4/12, 200e-4*200e-4/12,
1576 300e-4*300e-4/12, 300e-4*300e-4/12,
1577 300e-4*300e-4/12, 300e-4*300e-4/12}; // thickness^2/12
1578 //
1579 fTrack = NULL;
8102b2c9 1580 Int_t idx[20] = {0};
1581 Short_t lrID[20] = {0};
6be22b3f 1582 Int_t npts=atp->GetNPoints();
ef24eb3b 1583 if (npts<fMinNPtsPerTrack) return NULL;
6be22b3f 1584 TGeoHMatrix hcov;
1585 //
1586 /// checks if AliTrackPoints belong to defined modules
1587 Int_t ngoodpts=0;
1588 Int_t intidx[20];
1589 for (int j=0; j<npts; j++) {
8fd71c0a 1590 intidx[j] = GetRequestedModID(atp->GetVolumeID()[j]);
6be22b3f 1591 if (intidx[j]<0) continue;
1592 ngoodpts++;
1593 Float_t xx=atp->GetX()[j];
1594 Float_t yy=atp->GetY()[j];
1595 Float_t r=xx*xx + yy*yy;
1596 int lay;
1597 for (lay=0;lay<6;lay++) if (r<kRad2L[lay]) break;
1598 if (lay>5) continue;
1599 lrID[j] = lay;
1600 }
1601 //
1602 AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
1603
1604 // pepo270809
1605 Int_t nextra=0;
1606 // extra clusters selection mode
1607 if (fExtraClustersMode) {
1608 // 1 = keep one cluster, remove randomly the extra
1609 // 2 = keep one cluster, remove the internal one
1610 // 10 = keep tracks only if at least one extra is present
1611
1612 int iextra1[20],iextra2[20],layovl[20];
1613 // extra clusters mapping
1614 for (Int_t ipt=0; ipt<npts; ipt++) {
1615 if (intidx[ipt]<0) continue; // looks only defined modules...
1616 float p1x=atp->GetX()[ipt];
1617 float p1y=atp->GetY()[ipt];
1618 float p1z=atp->GetZ()[ipt];
1619 int lay1=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ipt]));
1620 float r1 = p1x*p1x + p1y*p1y;
1621 UShort_t volid1=atp->GetVolumeID()[ipt];
1622
1623 for (int ik=ipt+1; ik<npts; ik++) {
1624 if (intidx[ik]<0) continue;
1625 // compare point ipt with next ones
1626 int lay2=int(AliGeomManager::VolUIDToLayer(atp->GetVolumeID()[ik]));
1627 // check if same layer
1628 if (lay2 != lay1) continue;
1629 UShort_t volid2=atp->GetVolumeID()[ik];
1630 // check if different module
1631 if (volid1 == volid2) continue;
1632
1633 float p2x=atp->GetX()[ik];
1634 float p2y=atp->GetY()[ik];
1635 float p2z=atp->GetZ()[ik];
1636 float r2 = p2x*p2x + p2y*p2y;
1637 float dr= (p1x-p2x)*(p1x-p2x) + (p1y-p2y)*(p1y-p2y) + (p1z-p2z)*(p1z-p2z);
1638
1639 // looks for pairs with dr<1 cm, same layer but different module
1640 if (dr<1.0) {
1641 // extra1 is the one with smaller radius in rphi plane
1642 if (r1<r2) {
1643 iextra1[nextra]=ipt;
1644 iextra2[nextra]=ik;
1645 }
1646 else {
1647 iextra1[nextra]=ik;
1648 iextra2[nextra]=ipt;
1649 }
1650 layovl[nextra]=lay1;
1651 nextra++;
1652 }
1653 }
1654 } // end overlaps mapping
1655
1656 // mode=1: keep only one clusters and remove the other randomly
1657 if (fExtraClustersMode==1 && nextra) {
1658 for (int ie=0; ie<nextra; ie++) {
1659 if (gRandom->Rndm()<0.5)
1660 intidx[iextra1[ie]]=-1;
1661 else
1662 intidx[iextra2[ie]]=-1;
1663 }
1664 }
1665
1666 // mode=2: keep only one clusters and remove the other...
1667 if (fExtraClustersMode==2 && nextra) {
1668 for (int ie=0; ie<nextra; ie++) {
1669 if (layovl[ie]==1) intidx[iextra2[ie]]=-1;
1670 else if (layovl[ie]==2) intidx[iextra1[ie]]=-1;
1671 else intidx[iextra1[ie]]=-1;
1672 }
1673 }
1674
1675 // mode=10: reject track if no overlaps are present
1676 if (fExtraClustersMode==10 && nextra==0) {
1677 AliInfo("Track with no extra clusters: rejected!");
1678 return NULL;
1679 }
1680
1681 // recalculate ngoodpts
1682 ngoodpts=0;
1683 for (int i=0; i<npts; i++) {
1684 if (intidx[i]>=0) ngoodpts++;
1685 }
1686 }
1687 // endpepo270809
1688
1689 // reject track if not enough points are left
1690 if (ngoodpts<fMinNPtsPerTrack) {
ef24eb3b 1691 AliDebug(2,"Track with not enough points!");
6be22b3f 1692 return NULL;
1693 }
1694 // >> RS
1695 AliTrackPoint p;
1696 // check points in specific places
1d06ac63 1697 if (fRequirePoints[fDataType]) {
6be22b3f 1698 Int_t nlayup[6],nlaydown[6],nlay[6];
1699 Int_t ndetup[3],ndetdown[3],ndet[3];
1700 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1701 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1702
1703 for (int i=0; i<npts; i++) {
1704 // skip not defined points
1705 if (intidx[i]<0) continue;
1706 //
1707 Float_t yy=atp->GetY()[i];
1708 int lay = lrID[i];
1709 int det=lay/2;
1710 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
1711
1712 if (yy>=0.0) { // UP point
1713 nlayup[lay]++;
1714 nlay[lay]++;
1715 ndetup[det]++;
1716 ndet[det]++;
1717 }
1718 else {
1719 nlaydown[lay]++;
1720 nlay[lay]++;
1721 ndetdown[det]++;
1722 ndet[det]++;
1723 }
1724 }
1725 //
1726 // checks minimum values
1727 Bool_t isok=kTRUE;
1728 for (Int_t j=0; j<6; j++) {
1d06ac63 1729 if (nlayup[j]<fNReqLayUp[fDataType][j]) isok=kFALSE;
1730 if (nlaydown[j]<fNReqLayDown[fDataType][j]) isok=kFALSE;
1731 if (nlay[j]<fNReqLay[fDataType][j]) isok=kFALSE;
6be22b3f 1732 }
1733 for (Int_t j=0; j<3; j++) {
1d06ac63 1734 if (ndetup[j]<fNReqDetUp[fDataType][j]) isok=kFALSE;
1735 if (ndetdown[j]<fNReqDetDown[fDataType][j]) isok=kFALSE;
1736 if (ndet[j]<fNReqDet[fDataType][j]) isok=kFALSE;
6be22b3f 1737 }
1738 if (!isok) {
1739 AliDebug(2,Form("Track does not meet all location point requirements!"));
1740 return NULL;
1741 }
1742 }
1743 // build a new track with (sorted) (prealigned) good points
1744 // pepo200709
1745 //fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
8102b2c9 1746 Int_t addVertex = IsTypeCollision()&&((fUseDiamond&&(fCheckDiamondPoint!=kDiamondIgnore))||(fUseVertex&&fVertexSet)) ? 1 : 0;
ef24eb3b 1747 fTrack = (AliTrackPointArray*)fTrackBuff[ngoodpts + addVertex ];
6be22b3f 1748 if (!fTrack) {
ef24eb3b 1749 fTrack = new AliTrackPointArray(ngoodpts + addVertex);
6be22b3f 1750 // fTrackBuff.AddAtAndExpand(fTrack,ngoodpts-fMinNPtsPerTrack);
ef24eb3b 1751 fTrackBuff.AddAtAndExpand(fTrack,ngoodpts + addVertex);
6be22b3f 1752 }
1753 // fTrack = new AliTrackPointArray(ngoodpts);
1754 // endpepo200709
1755 //
1756 //
1757 for (int i=0; i<npts; i++) idx[i]=i;
1758 // sort track if required
8102b2c9 1759 if (IsTypeCosmics()) TMath::Sort(npts,atp->GetY(),idx); // sort descending...
6be22b3f 1760 //
1761 Int_t npto=0;
1762 if (fClusLoc.GetSize()<3*npts) fClusLoc.Set(3*npts);
1763 if (fClusGlo.GetSize()<3*npts) fClusGlo.Set(3*npts);
1764 if (fClusSigLoc.GetSize()<3*npts) fClusSigLoc.Set(3*npts);
1765 //
1766 for (int i=0; i<npts; i++) {
1767 // skip not defined points
1768 if (intidx[idx[i]]<0) continue;
1769 atp->GetPoint(p,idx[i]);
1770 int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1771 //
1772 // prealign point if required
1773 // get matrix used to produce the digits
1774 AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1775 TGeoHMatrix *svOrigMatrix = GetSensorOrigMatrixSID(sid); //mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1776 // get back real local coordinate
ef24eb3b 1777 fMeasLoc = fClusLoc.GetArray() + npto*3;
1778 fMeasGlo = fClusGlo.GetArray() + npto*3;
1779 fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1780 fMeasGlo[0]=p.GetX();
1781 fMeasGlo[1]=p.GetY();
1782 fMeasGlo[2]=p.GetZ();
1783 AliDebug(3,Form("Global coordinates of measured point : X=%+f Y=%+f Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
1784 svOrigMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1785 AliDebug(3,Form("Local coordinates of measured point : X=%+f Y=%+f Z=%+f \n",fMeasLoc[0],fMeasLoc[1],fMeasLoc[2]));
1786 //
1787 if (p.GetDriftTime()>0) ProcessSDDPointInfo(&p,sid, npto); // for SDD points extract vdrift
1788 //
6be22b3f 1789 // update covariance matrix
1790 Double_t hcovel[9];
1791 hcovel[0]=double(p.GetCov()[0]);
1792 hcovel[1]=double(p.GetCov()[1]);
1793 hcovel[2]=double(p.GetCov()[2]);
1794 hcovel[3]=double(p.GetCov()[1]);
1795 hcovel[4]=double(p.GetCov()[3]);
1796 hcovel[5]=double(p.GetCov()[4]);
1797 hcovel[6]=double(p.GetCov()[2]);
1798 hcovel[7]=double(p.GetCov()[4]);
1799 hcovel[8]=double(p.GetCov()[5]);
1800 hcov.SetRotation(hcovel);
abd7ef79 1801 //
1802 if (AliLog::GetGlobalDebugLevel()>=2) {
1803 AliInfo("Original Global Cov Matrix");
1804 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovel[0],hcovel[1],hcovel[2],hcovel[4],hcovel[5],hcovel[8]);
1805 }
1806 //
6be22b3f 1807 // now rotate in local system
1808 hcov.Multiply(svOrigMatrix);
1809 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1810 // now hcov is LOCAL COVARIANCE MATRIX
1811 // apply sigma scaling
abd7ef79 1812 Double_t *hcovscl = hcov.GetRotationMatrix();
8102b2c9 1813 /*
1814 const float *cv = p.GetCov();
1815 printf("## %d %d %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e %+.3e\n",
1816 sid,p.GetClusterType(),
1817 fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],
1818 fMeasLoc[0],fMeasLoc[1],fMeasLoc[2],
1819 cv[0],cv[1],cv[2],cv[3],cv[4],cv[5],
1820 hcovscl[0],hcovscl[4],hcovscl[8]);
1821
1822 */
abd7ef79 1823 if (AliLog::GetGlobalDebugLevel()>=2) {
1824 AliInfo("Original Local Cov Matrix");
1825 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1826 }
1827 hcovscl[4] = fUseLocalYErr ? kSensSigY2[lrID[idx[i]]] : 1E-8; // error due to the sensor thickness
6be22b3f 1828 //
1829 for (int ir=3;ir--;) for (int ic=3;ic--;) {
1830 if (ir==ic) {
abd7ef79 1831 if ( IsZero(hcovscl[ir*3+ic],1e-8) ) hcovscl[ir*3+ic] = 1E-8;
6be22b3f 1832 else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
ef24eb3b 1833 fSigmaLoc[ir] = TMath::Sqrt(hcovscl[ir*3+ic]);
6be22b3f 1834 }
1835 else hcovscl[ir*3+ic] = 0;
1836 }
1837 //
abd7ef79 1838 if (AliLog::GetGlobalDebugLevel()>=2) {
1839 AliInfo("Modified Local Cov Matrix");
1840 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1841 }
6be22b3f 1842 //
1843 if (fBug==1) {
1844 // correzione bug LAYER 5 SSD temporanea..
1845 int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1846 if (ssdidx>=500 && ssdidx<1248) {
1847 int ladder=(ssdidx-500)%22;
1848 if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1849 if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1850 }
1851 }
1852 /// get (evenctually prealigned) matrix of sens. vol.
1853 TGeoHMatrix *svMatrix = GetSensorCurrMatrixSID(sid); //mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
1854 // modify global coordinates according with pre-aligment
ef24eb3b 1855 svMatrix->LocalToMaster(fMeasLoc,fMeasGlo);
6be22b3f 1856 // now rotate in local system
1857 hcov.Multiply(&svMatrix->Inverse());
1858 hcov.MultiplyLeft(svMatrix); // hcov is back in GLOBAL RF
1859 // cure once more
1860 for (int ir=3;ir--;) for (int ic=3;ic--;) if (IsZero(hcovscl[ir*3+ic])) hcovscl[ir*3+ic] = 0.;
1861 // printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1862 //
abd7ef79 1863 if (AliLog::GetGlobalDebugLevel()>=2) {
1864 AliInfo("Modified Global Cov Matrix");
1865 printf("%+.4e %+.4e %+.4e\n%+.4e %+.4e\n%+.4e\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[4],hcovscl[5],hcovscl[8]);
1866 }
1867 //
6be22b3f 1868 Float_t pcov[6];
1869 pcov[0]=hcovscl[0];
1870 pcov[1]=hcovscl[1];
1871 pcov[2]=hcovscl[2];
1872 pcov[3]=hcovscl[4];
1873 pcov[4]=hcovscl[5];
1874 pcov[5]=hcovscl[8];
8102b2c9 1875 //
1876 // make sure the matrix is positive definite
1877 {
1878 enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1879 if (pcov[kXX]*pcov[kYY]*0.999<pcov[kXY]*pcov[kXY]) pcov[kXY] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kYY]),pcov[kXY]);
1880 if (pcov[kXX]*pcov[kZZ]*0.999<pcov[kXZ]*pcov[kXZ]) pcov[kXZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kXX]*pcov[kZZ]),pcov[kXZ]);
1881 if (pcov[kYY]*pcov[kZZ]*0.999<pcov[kYZ]*pcov[kYZ]) pcov[kYZ] = 0.999*TMath::Sign((float)TMath::Sqrt(pcov[kYY]*pcov[kZZ]),pcov[kYZ]);
1882 }
6be22b3f 1883 //
ef24eb3b 1884 p.SetXYZ(fMeasGlo[0],fMeasGlo[1],fMeasGlo[2],pcov);
1885 // printf("New Gl coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]);
1886 AliDebug(3,Form("New global coordinates of measured point : X=%+f Y=%+f Z=%+f \n",fMeasGlo[0],fMeasGlo[1],fMeasGlo[2]));
6be22b3f 1887 fTrack->AddPoint(npto,&p);
ef24eb3b 1888 AliDebug(2,Form("Adding point[%d] = ( %+f , %+f , %+f ) volid = %d",npto,fTrack->GetX()[npto],
6be22b3f 1889 fTrack->GetY()[npto],fTrack->GetZ()[npto],fTrack->GetVolumeID()[npto] ));
1890 // printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY());
1891 npto++;
1892 }
1893 //
ef24eb3b 1894 fDiamondPointID = -1;
1895 if (addVertex) {
1896 fTrack->AddPoint(npto,&fDiamond);
1897 fMeasLoc = fClusLoc.GetArray() + npto*3;
1898 fMeasGlo = fClusGlo.GetArray() + npto*3;
1899 fSigmaLoc = fClusSigLoc.GetArray() + npto*3;
1900 fMeasLoc[0] = fMeasGlo[0] = fDiamond.GetX();
1901 fMeasLoc[1] = fMeasGlo[1] = fDiamond.GetY();
1902 fMeasLoc[2] = fMeasGlo[2] = fDiamond.GetZ();
8102b2c9 1903 fSigmaLoc[0] = TMath::Sqrt(fDiamond.GetCov()[0]);
1904 fSigmaLoc[1] = TMath::Sqrt(fDiamond.GetCov()[3]);
1905 fSigmaLoc[2] = TMath::Sqrt(fDiamond.GetCov()[5]);
ef24eb3b 1906 fDiamondPointID = npto++;
1907 }
1908 //
6be22b3f 1909 return fTrack;
1910}
1911
1912//________________________________________________________________________________________________________
1913AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp)
1914{
1915 /// sort alitrackpoints w.r.t. global Y direction
1916 AliTrackPointArray *atps=NULL;
1917 Int_t idx[20];
1918 Int_t npts=atp->GetNPoints();
1919 AliTrackPoint p;
1920 atps=new AliTrackPointArray(npts);
1921
1922 TMath::Sort(npts,atp->GetY(),idx);
1923
1924 for (int i=0; i<npts; i++) {
1925 atp->GetPoint(p,idx[i]);
1926 atps->AddPoint(i,&p);
ef24eb3b 1927 AliDebug(2,Form("Point[%d] = ( %+f , %+f , %+f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
6be22b3f 1928 }
1929 return atps;
1930}
1931
1932//________________________________________________________________________________________________________
1933Int_t AliITSAlignMille2::GetCurrentLayer() const
1934{
1935 // get current layer id
1936 if (!fGeoManager) {
1937 AliInfo("ITS geometry not initialized!");
1938 return -1;
1939 }
1940 return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1941}
1942
1943//________________________________________________________________________________________________________
1944Int_t AliITSAlignMille2::InitModuleParams()
1945{
1946 /// initialize geometry parameters for a given detector
1947 /// for current cluster (fCluster)
1948 /// fGlobalInitParam[] is set as:
1949 /// [tx,ty,tz,psi,theta,phi]
1950 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1951 /// *** At the moment: using Raffalele's angles definition ***
1952 ///
1953 /// return 0 if success
1954 /// If module is found but has no parameters to vary, return 1
1955
1956 if (!fGeoManager) {
1957 AliInfo("ITS geometry not initialized!");
1958 return -1;
1959 }
1960
1961 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1962
1963 // set the internal index (index in module list)
1964 UShort_t voluid=fCluster.GetVolumeID();
1965 fCurrentSensID = AliITSAlignMille2Module::GetIndexFromVolumeID(voluid);
1966 //
ef24eb3b 1967 if (fCurrentSensID==-1) { // this is a special "vertex" module
1968 fCurrentModule = GetMilleModuleByVID(voluid);
1969 fCurrentSensID = fCurrentModule->GetIndex();
1970
1971 }
1972 else {
1973 //
1974 // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1975 Int_t k=fNModules-1;
1976 fCurrentModule = 0;
1977 // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules
1978 while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
1979 if (k<0) return -3;
1980 }
6be22b3f 1981 //
1982 for (int i=AliITSAlignMille2Module::kMaxParTot;i--;) fModuleInitParam[i] = 0.0;
1983 //
1984 int clID = fCluster.GetUniqueID()-1;
1985 if (clID<0) { // external cluster
1986 fMeasGlo = &fExtClusterPar[0];
1987 fMeasLoc = &fExtClusterPar[3];
1988 fSigmaLoc = &fExtClusterPar[6];
1989 fExtClusterPar[0] = fCluster.GetX();
1990 fExtClusterPar[1] = fCluster.GetY();
1991 fExtClusterPar[2] = fCluster.GetZ();
1992 //
1993 TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1994 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1995 TGeoHMatrix hcov;
1996 Double_t hcovel[9];
1997 hcovel[0]=double(fCluster.GetCov()[0]);
1998 hcovel[1]=double(fCluster.GetCov()[1]);
1999 hcovel[2]=double(fCluster.GetCov()[2]);
2000 hcovel[3]=double(fCluster.GetCov()[1]);
2001 hcovel[4]=double(fCluster.GetCov()[3]);
2002 hcovel[5]=double(fCluster.GetCov()[4]);
2003 hcovel[6]=double(fCluster.GetCov()[2]);
2004 hcovel[7]=double(fCluster.GetCov()[4]);
2005 hcovel[8]=double(fCluster.GetCov()[5]);
2006 hcov.SetRotation(hcovel);
2007 // now rotate in local system
2008 hcov.Multiply(svMatrix);
2009 hcov.MultiplyLeft(&svMatrix->Inverse());
2010 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2011 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
2012 //
2013 }
2014 else {
2015 int offs = 3*clID;
2016 fMeasGlo = fClusGlo.GetArray() + offs;
2017 fMeasLoc = fClusLoc.GetArray() + offs;
2018 fSigmaLoc = fClusSigLoc.GetArray() + offs;
2019 }
2020 //
2021 // set minimum value for SigmaLoc to 10 micron
2022 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
2023 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
8102b2c9 2024 if (fCurrentSensID==kVtxSensID || fUseLocalYErr) if (fSigmaLoc[1]<0.0010) fSigmaLoc[1]=0.0010;
6be22b3f 2025 //
ef24eb3b 2026 AliDebug(2,Form("Local coordinates of measured point : X=%+f Y=%+f Z=%+f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
6be22b3f 2027 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
2028 //
2029 return 0;
2030}
2031
2032//________________________________________________________________________________________________________
2033void AliITSAlignMille2::Print(Option_t*) const
2034{
2035 // print current status
2036 printf("*** AliMillepede for ITS ***\n");
2037 printf(" Number of defined super modules: %d\n",fNModules);
2038 printf(" Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
2039 //
2040 if (fGeoManager)
2041 printf(" geometry loaded from %s\n",fGeometryPath.Data());
2042 else
2043 printf(" geometry not loaded\n");
2044 //
2045 if (fUsePreAlignment)
2046 printf(" using prealignment from %s \n",fPreDeltaPath.Data());
2047 else
2048 printf(" prealignment not used\n");
2049 //
2050 //
2051 if (fBOn)
ef24eb3b 2052 printf(" B Field set to %+f T - using helices\n",fBField);
6be22b3f 2053 else
2054 printf(" B Field OFF - using straight lines \n");
2055 //
2056 if (fTPAFitter)
2057 printf(" Using AliITSTPArrayFit class for track fitting\n");
2058 else
2059 printf(" Using StraightLine/Riemann fitter for track fitting\n");
2060 //
2061 printf("Using local Y error due to the sensor thickness: %s\n",(fUseLocalYErr && fTPAFitter) ? "ON":"OFF");
2062 //
1d06ac63 2063 for (int itp=0;itp<kNDataType;itp++) {
2064 if (fRequirePoints[itp]) printf(" Required points in %s tracks:\n",itp==kCosmics? "cosmics" : "collisions");
2065 for (Int_t i=0; i<6; i++) {
2066 if (fNReqLayUp[itp][i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[itp][i]);
2067 if (fNReqLayDown[itp][i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[itp][i]);
2068 if (fNReqLay[itp][i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[itp][i]);
2069 }
2070 for (Int_t i=0; i<3; i++) {
2071 if (fNReqDetUp[itp][i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[itp][i]);
2072 if (fNReqDetDown[itp][i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[itp][i]);
2073 if (fNReqDet[itp][i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[itp][i]);
2074 }
6be22b3f 2075 }
ef24eb3b 2076 printf(" SDD VDrift correction : %s",fIsSDDVDriftMult ? "Mult":"Add");
2077 printf(" Weight acc. to pT in power : %f",fWeightPt);
6be22b3f 2078 //
2079 printf("\n Millepede configuration parameters:\n");
ef24eb3b 2080 printf(" init factor for chi2 cut : %.4f\n",fStartFac);
2081 printf(" final factor for chi2 cut : %.4f\n",fFinalFac);
6be22b3f 2082 printf(" first iteration cut value : %.4f\n",fResCutInitial);
2083 printf(" other iterations cut value : %.4f\n",fResCut);
2084 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
2085 printf(" def.scaling for local sigmas : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
2086 printf(" min.tracks per module : %d\n",fMinPntPerSens);
2087 //
2088 printf("List of defined modules:\n");
2089 printf(" intidx\tindex\tvoluid\tname\n");
2090 for (int i=0; i<fNModules; i++) {
2091 AliITSAlignMille2Module* md = GetMilleModule(i);
2092 printf(" %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
2093 }
2094}
2095
2096//________________________________________________________________________________________________________
2097AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
2098{
2099 // return pointer to a defined supermodule
2100 // return NULL if error
2101 Int_t i=IsVIDDefined(voluid);
2102 if (i<0) return NULL;
2103 return GetMilleModule(i);
2104}
2105
2106//________________________________________________________________________________________________________
2107AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
2108{
2109 // return pointer to a defined supermodule
2110 // return NULL if error
2111 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2112 if (vid>0) return GetMilleModuleByVID(vid);
2113 else { // this is not alignable module, need to look within defined supermodules
2114 int i = IsSymDefined(symname);
2115 if (i>=0) return GetMilleModule(i);
2116 }
2117 return 0;
2118}
2119
2120//________________________________________________________________________________________________________
2121AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
2122{
2123 // return pointer to a defined/contained supermodule
2124 // return NULL otherwise
2125 int i = IsSymContained(symname);
2126 return i<0 ? 0 : GetMilleModule(i);
2127}
2128
2129//________________________________________________________________________________________________________
2130AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
2131{
2132 // get delta from prealignment for given volume
2133 if (!fPrealignment) return 0;
2134 for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2135 AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
2136 if (!strcmp(preob->GetSymName(),symname)) return preob;
2137 }
2138 return 0;
2139}
2140
2141//________________________________________________________________________________________________________
2142AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
2143{
2144 // get delta with respect to which the constraint is declared
2145 if (!fConstrRef) return 0;
2146 for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
2147 AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
2148 if (!strcmp(preob->GetSymName(),symname)) return preob;
2149 }
2150 return 0;
2151}
2152
2153//________________________________________________________________________________________________________
2154Bool_t AliITSAlignMille2::InitRiemanFit()
2155{
2156 // Initialize Riemann Fitter for current track
2157 // return kFALSE if error
2158
2159 if (!fBOn) return kFALSE;
2160
2161 Int_t npts=0;
2162 AliTrackPoint ap;
2163 npts = fTrack->GetNPoints();
2164 AliDebug(3,Form("Fitting track with %d points",npts));
2165 if (!fRieman) fRieman = new AliTrackFitterRieman();
2166 fRieman->Reset();
2167 fRieman->SetTrackPointArray(fTrack);
2168
2169 TArrayI ai(npts);
2170 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
2171
2172 // fit track with 5 params in his own tracking-rotated reference system
2173 // xc = -p[1]/p[0];
2174 // yc = 1/p[0];
2175 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
2176 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
2177 return kFALSE;
2178 }
2179
2180 for (int i=0; i<5; i++)
2181 fLocalInitParam[i] = fRieman->GetParam()[i];
2182
2183 return kTRUE;
2184}
2185
2186//________________________________________________________________________________________________________
2187void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int flag)
2188{
2189 // local function for minuit
2190 const double kTiny = 1.e-14;
2191 chi2 = 0;
2192 static AliTrackPoint pnt;
2193 static Bool_t fullErr2D;
2194 //
2195 if (flag==1) fullErr2D = kFALSE;//kTRUE;
8fd71c0a 2196 // fullErr2D = kTRUE;
6be22b3f 2197 enum {kAX,kAZ,kBX,kBZ};
2198 enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
2199 //
2200 AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
2201 AliTrackPointArray* track = alig->GetCurrentTrack();
2202 //
2203 int npts = track->GetNPoints();
2204 for (int ip=0;ip<npts;ip++) {
2205 track->GetPoint(pnt,ip);
2206 const float *cov = pnt.GetCov();
2207 double y = pnt.GetY();
2208 double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
2209 double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
2210 double xxe = cov[kXX];
2211 double zze = cov[kZZ];
2212 double xze = cov[kXZ];
2213 //
2214 if (fullErr2D) {
2215 xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
2216 zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
2217 xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
2218 }
2219 //
2220 double det = xxe*zze - xze*xze;
2221 if (det<kTiny) {
2222 printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
2223 "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
2224 xxe = cov[kXX];
2225 zze = cov[kZZ];
2226 xze = cov[kXZ];
2227 fullErr2D = kFALSE;
2228 }
2229 double xxeI = zze/det;
2230 double zzeI = xxe/det;
2231 double xzeI =-xze/det;
2232 //
2233 chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
2234 //
2235 // printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI, chi2);
2236 }
2237 //
2238}
2239
2240//________________________________________________________________________________________________________
2241void AliITSAlignMille2::InitTrackParams(int meth)
2242{
2243 /// initialize local parameters with different methods
2244 /// for current track (fTrack)
2245 Int_t npts=0;
2246 AliTrackPoint ap;
2247 double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
2248 // simple linear interpolation
2249 // get local starting parameters (to be substituted by ESD track parms)
2250 // local parms (fLocalInitParam[]) are:
2251 // [0] = global x coord. of straight line intersection at y=0 plane
2252 // [1] = global z coord. of straight line intersection at y=0 plane
2253 // [2] = px/py
2254 // [3] = pz/py
2255 // test #1: linear fit in x(y) and z(y)
2256 npts = fTrack->GetNPoints();
2257 AliDebug(3,Form("*** initializing track with %d points ***",npts));
2258 for (int i=npts;i--;) {
2259 sY += fTrack->GetY()[i];
2260 sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
2261 sX += fTrack->GetX()[i];
2262 sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
2263 sZ += fTrack->GetZ()[i];
2264 sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
2265 }
2266 det = sYY*npts-sY*sY;
2267 if (IsZero(det)) det = 1E-16;
2268 fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
2269 fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
2270 //
2271 fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
2272 fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
2273 // pepo200709
2274 fLocalInitParam[4] = 0.0;
2275 // endpepo200709
2276
ef24eb3b 2277 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %+f ugx = %+f\n",fLocalInitParam[0],fLocalInitParam[2]));
6be22b3f 2278 //
2279 if (meth==1) return;
2280 //
2281 // perform full fit accounting for cov.matrix
2282 static TVirtualFitter *minuit = 0;
2283 static Double_t step[5] = {1E-3,1E-3,1E-4,1E-4,1E-5};
2284 static Double_t arglist[10];
2285 //
2286 if (!minuit) {
2287 minuit = TVirtualFitter::Fitter(0,4);
2288 minuit->SetFCN(trackFit2D);
2289 arglist[0] = 1;
2290 minuit->ExecuteCommand("SET ERR",arglist, 1);
2291 //
2292 arglist[0] = -1;
2293 minuit->ExecuteCommand("SET PRINT",arglist,1);
2294 //
2295 }
2296 //
2297 minuit->SetParameter(0, "ax", fLocalInitParam[0], step[0], 0,0);
2298 minuit->SetParameter(1, "az", fLocalInitParam[1], step[1], 0,0);
2299 minuit->SetParameter(2, "bx", fLocalInitParam[2], step[2], 0,0);
2300 minuit->SetParameter(3, "bz", fLocalInitParam[3], step[3], 0,0);
2301 //
2302 arglist[0] = 1000; // number of function calls
2303 arglist[1] = 0.001; // tolerance
2304 minuit->ExecuteCommand("MIGRAD",arglist,2);
2305 //
2306 for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
2307 for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
2308 /*
2309 double amin,edm,errdef;
2310 int nvpar,nparx;
2311 minuit->GetStats(amin,edm,errdef,nvpar,nparx);
2312 amin /= (2*npts - 4);
2313 printf("Mchi2: %+e\n",amin);
2314 */
2315 //
2316}
2317
2318//________________________________________________________________________________________________________
2319Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
2320{
2321 // checks if supermodule with this symname is defined and return the internal index
2322 // return -1 if not.
2323 for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
2324 return -1;
2325}
2326
2327//________________________________________________________________________________________________________
2328Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
2329{
2330 // checks if module with this symname is defined and return the internal index
2331 // return -1 if not.
2332 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
2333 if (vid>0) return IsVIDContained(vid);
2334 // only sensors have real vid, but maybe we have a supermodule with fake vid?
2335 // IMPORTANT: always start from the end to start from the sensors
2336 return IsSymDefined(symname);
2337}
2338
2339//________________________________________________________________________________________________________
2340Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
2341{
2342 // checks if supermodule 'voluid' is defined and return the internal index
2343 // return -1 if not.
2344 for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
2345 return -1;
2346}
2347
2348//________________________________________________________________________________________________________
2349Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
2350{
2351 // checks if the sensitive module 'voluid' is contained inside a supermodule
2352 // and return the internal index of the last identified supermodule
2353 // return -1 if error
2354 // IMPORTANT: always start from the end to start from the sensors
2355 if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2356 for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
2357 return -1;
2358}
2359
8fd71c0a 2360//________________________________________________________________________________________________________
2361Int_t AliITSAlignMille2::GetRequestedModID(UShort_t voluid) const
2362{
2363 // checks if the sensitive module 'voluid' is contained inside a supermodule
2364 // and return the internal index of the last identified supermodule
2365 // return -1 if error
2366 // IMPORTANT: always start from the end to start from the sensors
2367 if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
2368 int k;
2369 for (k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) break;
2370 if (k<0) return -1;
2371 AliITSAlignMille2Module* md = GetMilleModule(k);
2372 while (md && md->IsNotInConf()) md = md->GetParent();
3dc409f6 2373 if (md) return int(md->GetUniqueID());
2374 else return -1;
8fd71c0a 2375}
2376
6be22b3f 2377//________________________________________________________________________________________________________
2378Int_t AliITSAlignMille2::CheckCurrentTrack()
2379{
2380 /// checks if AliTrackPoints belongs to defined modules
2381 /// return number of good poins
2382 /// return 0 if not enough points
2383
2384 Int_t npts = fTrack->GetNPoints();
2385 Int_t ngoodpts=0;
2386 // debug points
2387 for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
2388 //
2389 if (ngoodpts<fMinNPtsPerTrack) return 0;
2390
2391 return ngoodpts;
2392}
2393
2394//________________________________________________________________________________________________________
ef24eb3b 2395Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track, Double_t wgh)
6be22b3f 2396{
2397 /// Process track; Loop over hits and set local equations
2398 /// here 'track' is a AliTrackPointArray
2399 /// return 0 if success;
ef24eb3b 2400 //
6be22b3f 2401 if (!fIsMilleInit) Init();
2402 //
2403 Int_t npts = track->GetNPoints();
2404 AliDebug(2,Form("*** Input track with %d points ***",npts));
2405
2406 // preprocessing of the input track: keep only points in defined volumes,
2407 // move points if prealignment is set, sort by Yglo if required
ef24eb3b 2408 fTrackWeight = wgh;
6be22b3f 2409 fTrack=PrepareTrack(track);
1d06ac63 2410 if (!fTrack) {
2411 RemoveHelixFitConstraint();
8102b2c9 2412 RemoveVertexConstraint();
1d06ac63 2413 return -1;
2414 }
6be22b3f 2415 npts = fTrack->GetNPoints();
2416 if (npts>kMaxPoints) {
2417 AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
2418 }
2419 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
2420 //
ef24eb3b 2421 npts = FitTrack();
2422 if (npts<0) return npts;
2423 //
2424 // printf("Params: "); for (int i=0;i<fNLocal;i++) printf("%+.2e ",fLocalInitParam[i]); printf("\n");//RRR
2425 Int_t nloceq=0;
2426 Int_t ngloeq=0;
2427 static Mille2Data md[kMaxPoints];
2428 //
2429 for (Int_t ipt=0; ipt<npts; ipt++) {
2430 fTrack->GetPoint(fCluster,ipt);
2431 fCluster.SetUniqueID(ipt+1);
2432 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
2433
2434 // set geometry parameters for the the current module
2435 if (InitModuleParams()) continue;
2436 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n",
2437 track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
2438 fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
2439 AliDebug(2,Form(" Preprocessed Point = ( %+f , %+f , %+f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
2440 int res = fTPAFitter ? AddLocalEquationTPA(md[nloceq]) : AddLocalEquation(md[nloceq]);
2441 if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
2442 else if (res==0) nloceq++;
2443 else {nloceq++; ngloeq++;}
2444 } // end loop over points
2445 //
2446 fTrack=NULL;
2447 // not enough good points?
2448 if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
2449 //
2450 // finally send local equations to millepede
2451 SetLocalEquations(md,nloceq);
2452 fMillepede->SaveRecordData(); // RRR
8693b6b7 2453 fCurvFitWasConstrained = kFALSE; // restore default
ef24eb3b 2454 //
2455 return 0;
2456}
2457
2458//________________________________________________________________________________________________________
2459Int_t AliITSAlignMille2::FitTrack()
2460{
2461 // Fit the track with selected constraints
2462 //
2463 const Double_t kfDiamondTolerance = 0.1; //diamond tolerance on top of the MS error
2464 if (!fTrack) return -1;
2465 int npts = fTrack->GetNPoints();
2466 //
6be22b3f 2467 if (fTPAFitter) { // use dediacted fitter
2468 //
ef24eb3b 2469 // if the diamond point is attached, for the moment don't include it in the fit
8102b2c9 2470 fTPAFitter->AttachPoints(fTrack,0, npts-1);
1d06ac63 2471 fTPAFitter->SetBz(fBField);
2472 fTPAFitter->SetTypeCosmics(IsTypeCosmics());
ef24eb3b 2473 if (fIniTrackParamsMeth==1) fTPAFitter->SetIgnoreCov();
6be22b3f 2474 //
ef24eb3b 2475 double chi2;
2476 double chi2f = 0;
2477 double dca2err;
2478 double dca2 = 0.;
2479 Bool_t fitIsDone = kFALSE;
8102b2c9 2480 if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) { // diamond constraint was added, check if the track looks like prompt
2481 fTPAFitter->SetFirstLast(0,fDiamondPointID-1);
2482 if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2483 //
ef24eb3b 2484 chi2f = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2485 if ( chi2f<0 || (chi2f>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations()) ) { //RRR
2486 AliInfo(Form("Track fit failed on checking if it is prompt! skipping this track... Chi2:%+e",chi2f));
2487 fTPAFitter->Reset();
2488 // fTrack = NULL;
2489 return -5;
2490 }
2491 double xyzRes[3];
2492 fTPAFitter->GetResiduals(xyzRes,&fDiamondI,kTRUE);
2493 dca2 = xyzRes[0]*xyzRes[0] + xyzRes[1]*xyzRes[1];
2494 double pT = IsFieldON() ? fTPAFitter->GetPt() : 0.45;
2495 if (pT<0.1) pT = 0.1;
2496 dca2err = kfDiamondTolerance + 0.02/pT;
2497 if (dca2>dca2err*dca2err) { // this is secondary
2498 int* clst = (int*) fTrack->GetClusterType();
2499 clst[fDiamondPointID] = -1;;
2500 fDiamondPointID = -1;
2501 fitIsDone = kTRUE;
2502 npts--;
2503 }
2504 else fTPAFitter->SetFirstLast(0,fDiamondPointID); // fit with diamond
2505 }
2506 // fTPAFitter->SetParAxis(1);
8102b2c9 2507 if (!fitIsDone) {
2508 if (IsCovIScaleTouched()) for (int i=npts;i--;) fTPAFitter->SetCovIScale(i,GetCovIScale(i));
2509 chi2 = fTPAFitter->Fit(fConstrCharge,fConstrPT,fConstrPTErr);
2510 }
ef24eb3b 2511 //
2512 RemoveHelixFitConstraint(); // suppress eventual constraints to not affect fit of the next track
8102b2c9 2513 RemoveVertexConstraint(); // same ...
6be22b3f 2514 //
ef24eb3b 2515 if ( !fitIsDone && (chi2<0 || (chi2>fNStdDev*fStartFac && fTPAFitter->GetNIterations()>=fTPAFitter->GetMaxIterations())) ) { //RRR
2516 AliInfo(Form("Track fit failed! skipping this track... Chi2:%+e",chi2));
8102b2c9 2517 if (fUseDiamond && fDiamondPointID>0 && fCheckDiamondPoint==kDiamondCheckIfPrompt) AliInfo(Form("VertexFree fit gave Chi2:%+e with residual %+e",chi2f,TMath::Sqrt(dca2)));
ef24eb3b 2518 /*
2519 fTrack->Print("");
2520 fTPAFitter->FitHelixCrude();
2521 fTPAFitter->SetFitDone();
2522 fTPAFitter->Print();
2523 */
6be22b3f 2524 fTPAFitter->Reset();
ef24eb3b 2525 // fTrack = NULL;
6be22b3f 2526 return -5;
2527 }
ef24eb3b 2528 fNLocal = fTPAFitter->IsFieldON() ? 5:4; // Attention: the fitter might have decided to work in line mode
2529 npts = fTPAFitter->GetLast() - fTPAFitter->GetFirst() + 1; // actual number of points
6be22b3f 2530 /*
ef24eb3b 2531 double *pr = fTPAFitter->GetParams();
2532 printf("FtPar: %+.5e %+.5e %+.5e %+.5e | chi2:%.3e\n",pr[2],pr[0],pr[3],pr[1],chi2); // RRR
6be22b3f 2533 */
2534 }
2535 else {
2536 //
2537 if (!fBOn) { // straight lines
2538 // set local starting parameters (to be substituted by ESD track parms)
2539 // local parms (fLocalInitParam[]) are:
2540 // [0] = global x coord. of straight line intersection at y=0 plane
2541 // [1] = global z coord. of straight line intersection at y=0 plane
2542 // [2] = px/py
2543 // [3] = pz/py
ef24eb3b 2544 InitTrackParams(fIniTrackParamsMeth);
6be22b3f 2545 /*
2546 double *pr = fLocalInitParam;
2547 printf("FtPar: %+.5e %+.5e %+.5e %+.5e |\n",pr[0],pr[1],pr[2],pr[3]); // RRR
2548 */
2549 }
2550 else {
2551 // local parms (fLocalInitParam[]) are the Riemann Fitter params
2552 if (!InitRiemanFit()) {
2553 AliInfo("Riemann fit failed! skipping this track...");
2554 fTrack=NULL;
2555 return -5;
2556 }
2557 }
2558 }
ef24eb3b 2559 return npts;
6be22b3f 2560 //
6be22b3f 2561}
2562
2563//________________________________________________________________________________________________________
b80c197e 2564Int_t AliITSAlignMille2::CalcIntersectionPoint(const Double_t *lpar, const Double_t *gpar)
6be22b3f 2565{
2566 /// calculate track intersection point in local coordinates
2567 /// according with a given set of parameters (local(4) and global(6))
2568 /// and fill fPintLoc/Glo
2569 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
2570 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
2571 /// return 0 if success
2572
2573 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]));
2574 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]));
2575
2576
2577 // prepare the TGeoHMatrix
2578 TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
2579 !fUseGlobalDelta);
2580 if (!tempHMat) return -1;
2581
2582 Double_t v0g[3]; // vector with straight line direction in global coord.
2583 Double_t p0g[3]; // point of the straight line (glo)
2584
2585 if (fBOn) { // B FIELD!
2586 AliTrackPoint prf;
2587 for (int ip=0; ip<5; ip++)
2588 fRieman->SetParam(ip,lpar[ip]);
2589
2590 if (!fRieman->GetPCA(fCluster,prf)) {
2591 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
2592 return -3;
2593 }
2594 // now determine straight line passing tangent to fit curve at prf
2595 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
2596 // mo' P1=(X,Y,Z)_glo_prf
2597 // => (x,y,Z)_trk_prf ruotando di alpha...
2598 Double_t alpha=fRieman->GetAlpha();
2599 Double_t x1g = prf.GetX();
2600 Double_t y1g = prf.GetY();
2601 Double_t z1g = prf.GetZ();
2602 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
2603 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
2604 Double_t z1t = z1g;
2605
2606 Double_t x2t = x1t+1.0;
2607 Double_t y2t = y1t+fRieman->GetDYat(x1t);
2608 Double_t z2t = z1t+fRieman->GetDZat(x1t);
2609 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
2610 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
2611 Double_t z2g = z2t;
2612
ef24eb3b 2613 AliDebug(3,Form("Riemann frame: fAlpha = %+f = %+f ",alpha,alpha*180./TMath::Pi()));
2614 AliDebug(3,Form(" prf_glo=( %+f , %+f , %+f ) prf_rf=( %+f , %+f , %+f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
2615 AliDebug(3,Form(" mov_glo=( %+f , %+f , %+f ) rf=( %+f , %+f , %+f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
6be22b3f 2616
2617 if (TMath::Abs(y2g-y1g)<1e-15) {
2618 AliInfo("DeltaY=0! Cannot proceed...");
2619 return -1;
2620 }
2621 // ugx,1,ugz
2622 v0g[0] = (x2g-x1g)/(y2g-y1g);
2623 v0g[1]=1.0;
2624 v0g[2] = (z2g-z1g)/(y2g-y1g);
2625
2626 // point: just keep prf
2627 p0g[0]=x1g;
2628 p0g[1]=y1g;
2629 p0g[2]=z1g;
2630 }
2631 else { // staight line
2632 // vector of initial straight line direction in glob. coord
2633 v0g[0]=lpar[2];
2634 v0g[1]=1.0;
2635 v0g[2]=lpar[3];
2636
2637 // intercept in yg=0 plane in glob coord
2638 p0g[0]=lpar[0];
2639 p0g[1]=0.0;
2640 p0g[2]=lpar[1];
2641 }
ef24eb3b 2642 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]));
6be22b3f 2643
2644 // same in local coord.
2645 Double_t p0l[3],v0l[3];
2646 tempHMat->MasterToLocalVect(v0g,v0l);
2647 tempHMat->MasterToLocal(p0g,p0l);
2648
2649 if (TMath::Abs(v0l[1])<1e-15) {
2650 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
2651 return -1;
2652 }
2653
2654 // local intersection point
2655 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
2656 fPintLoc[1] = 0;
2657 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
2658
2659 // global intersection point
2660 tempHMat->LocalToMaster(fPintLoc,fPintGlo);
ef24eb3b 2661 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]));
6be22b3f 2662
2663 return 0;
2664}
2665
2666//________________________________________________________________________________________________________
2667Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
2668{
2669 /// calculate numerically (ROOT's style) the derivatives for
2670 /// local X intersection and local Z intersection
2671 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
2672 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
2673 /// return 0 if success
2674
2675 // copy initial parameters
2676 Double_t lpar[kNLocal];
2677 Double_t gpar[kNParCh];
2678 Double_t *derivative;
2679 for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
2680 for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
2681
2682 // trial with fixed dpar...
2683 Double_t dpar = 0.;
2684
2685 if (islpar) { // track parameters
2686 //dpar=fLocalInitParam[paridx]*0.001;
2687 // set minimum dpar
2688 derivative = fDerivativeLoc[paridx];
2689 if (!fBOn) {
2690 if (paridx<3) dpar=1.0e-4; // translations
2691 else dpar=1.0e-6; // direction
2692 }
2693 else { // B Field
2694 // pepo: proviamo con 1/1000, poi evenctually 1/100...
2695 Double_t dfrac=0.01;
2696 switch(paridx) {
2697 case 0:
2698 // RMS cosmics: 1e-4
2699 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2700 break;
2701 case 1:
2702 // RMS cosmics: 0.2
2703 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2704 break;
2705 case 2:
2706 // RMS cosmics: 9
2707 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2708 break;
2709 case 3:
2710 // RMS cosmics: 7
2711 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2712 break;
2713 case 4:
2714 // RMS cosmics: 0.3
2715 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
2716 break;
2717 }
2718 }
2719 }
2720 else { // alignment global parameters
2721 derivative = fDerivativeGlo[paridx];
2722 //dpar=fModuleInitParam[paridx]*0.001;
2723 if (paridx<3) dpar=1.0e-4; // translations
2724 else dpar=1.0e-2; // angles
2725 }
2726
2727 AliDebug(3,Form("+++ using dpar=%g",dpar));
2728
2729 // calculate derivative ROOT's like:
2730 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2731 Double_t pintl1[3]; // f(x-h)
2732 Double_t pintl2[3]; // f(x-h/2)
2733 Double_t pintl3[3]; // f(x+h/2)
2734 Double_t pintl4[3]; // f(x+h)
2735
2736 // first values
2737 if (islpar) lpar[paridx] -= dpar;
2738 else gpar[paridx] -= dpar;
2739 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2740 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2741
2742 // second values
2743 if (islpar) lpar[paridx] += dpar/2;
2744 else gpar[paridx] += dpar/2;
2745 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2746 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2747
2748 // third values
2749 if (islpar) lpar[paridx] += dpar;
2750 else gpar[paridx] += dpar;
2751 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2752 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2753
2754 // fourth values
2755 if (islpar) lpar[paridx] += dpar/2;
2756 else gpar[paridx] += dpar/2;
2757 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2758 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2759
2760 Double_t h2 = 1./(2.*dpar);
2761 Double_t d0 = pintl4[0]-pintl1[0];
2762 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
2763 derivative[0] = h2*(4*d2 - d0)/3.;
2764 if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
2765
2766 d0 = pintl4[2]-pintl1[2];
2767 d2 = 2.*(pintl3[2]-pintl2[2]);
2768 derivative[2] = h2*(4*d2 - d0)/3.;
2769 if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
2770
2771 AliDebug(3,Form("\n+++ derivatives +++ \n"));
2772 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2773 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
2774
2775 return 0;
2776}
2777
2778//________________________________________________________________________________________________________
2779Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2780{
2781 /// Define local equation for current cluster in X and Z coor.
2782 /// and store them to memory
2783 /// return -1 in case of failure to build some equation
2784 /// 0 if no free global parameters were found but local eq is built
2785 /// 1 if both local and global eqs are built
2786 //
2787 // store first intersection point
2788 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
2789 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2790
ef24eb3b 2791 AliDebug(2,Form("Intersect. point: L( %+f , %+f , %+f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
6be22b3f 2792
2793 // calculate local derivatives numerically
2794 Bool_t zeroX = kTRUE;
2795 Bool_t zeroZ = kTRUE;
2796 //
2797 for (Int_t i=0; i<fNLocal; i++) {
2798 if (CalcDerivatives(i,kTRUE)) return -1;
2799 m.fDerLoc[i][kX] = fDerivativeLoc[i][0];
2800 m.fDerLoc[i][kZ] = fDerivativeLoc[i][2];
2801 if (zeroX) zeroX = IsZero(fDerivativeLoc[i][0]);
2802 if (zeroZ) zeroZ = IsZero(fDerivativeLoc[i][2]);
2803 }
2804 // for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
2805 //
2806 if (zeroX) {AliInfo("Skipping: zero local X derivatives!"); return -1;}
2807 if (zeroZ) {AliInfo("Skipping: zero local Z derivatives!"); return -1;}
2808 //
8fd71c0a 2809 int status = 0;
6be22b3f 2810 int ifill = 0;
2811 //
2812 AliITSAlignMille2Module* endModule = fCurrentModule;
2813 //
2814 zeroX = zeroZ = kTRUE;
2815 Bool_t dfDone[kNParCh];
2816 for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2817 m.fNModFilled = 0;
2818 //
2819 // special block for SDD derivatives
2820 Double_t jacobian[kNParChGeom];
2821 Int_t nmodTested = 0;
2822 //
2823 do {
2824 if (fCurrentModule->GetNParFree()==0) continue;
2825 nmodTested++;
2826 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2827 //
2828 if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2829 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2830 if (!dfDone[i]) {
2831 if (CalcDerivatives(i,kFALSE)) return -1;
2832 else {
2833 dfDone[i] = kTRUE;
2834 if (zeroX) zeroX = IsZero(fDerivativeGlo[i][0]);
2835 if (zeroZ) zeroZ = IsZero(fDerivativeGlo[i][2]);
2836 }
2837 }
2838 //
2839 m.fDerGlo[ifill][kX] = fDerivativeGlo[i][0];
2840 m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][2];
2841 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2842 }
2843 //
2844 // specific for special sensors
ef24eb3b 2845 Int_t sddLR = -1;
6be22b3f 2846 if ( fCurrentModule->IsSDD() &&
ef24eb3b 2847 (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0 ||
2848 // fCurrentModule->GetParOffset(sddLR = fMeasLoc[kX]>0 ?
2849 fCurrentModule->GetParOffset(sddLR = GetVDriftSDD()>0 ?
2850 AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR)>=0)
2851 ) {
6be22b3f 2852 //
2853 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2854 // where V0 and T are the nominal drift velocity, time and time0
2855 // and the dT0 and dV are the corrections:
2856 // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2857 // dX/dV = dX/dxloc * dxloc/dV = dX/dxloc * (T-T0)
2858 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2859 //
ef24eb3b 2860 if (!dfDone[AliITSAlignMille2Module::kDOFT0] || !dfDone[sddLR]) {
6be22b3f 2861 //
2862 double dXdxlocsens=0., dZdxlocsens=0.;
2863 //
2864 // if the current module is the sensor itself and we work with local params, then
2865 // we can directly take dX/dxloc_sens dZ/dxloc_sens
2866 if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2867 if (!dfDone[AliITSAlignMille2Module::kDOFTX]) {
2868 CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2869 dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2870 }
2871 dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2872 dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2873 }
2874 //
2875 else { // need to perform some transformations
2876 // fetch the jacobian of the transformation from the sensors local frame to the frame
2877 // where the parameters are defined:
2878 // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2879 if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2880 AliITSAlignMille2Module::kDOFTX, jacobian);
2881 // Local: dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2882 else fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2883 AliITSAlignMille2Module::kDOFTX, jacobian);
2884 //
2885 for (int j=0;j<kNParChGeom;j++) {
2886 // need global derivative even if the j-th param is locked
2887 if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2888 dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2889 dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2890 }
2891 }
2892 //
2893 if (zeroX) zeroX = IsZero(dXdxlocsens);
2894 if (zeroZ) zeroZ = IsZero(dZdxlocsens);
2895 //
2896 double vdrift = GetVDriftSDD();
2897 double tdrift = GetTDriftSDD();
2898 //
2899 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2900 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2901 dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2902 //
ef24eb3b 2903 double mltCorr = fIsSDDVDriftMult ? TMath::Abs(vdrift) : 1;
2904 fDerivativeGlo[sddLR][0] = -dXdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2905 fDerivativeGlo[sddLR][2] = -dZdxlocsens*mltCorr*TMath::Sign(tdrift,vdrift);
2906 dfDone[sddLR] = kTRUE;
6be22b3f 2907 //
2908 }
2909 //
2910 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2911 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2912 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2913 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2914 }
2915 //
ef24eb3b 2916 if (fCurrentModule->GetParOffset(sddLR)>=0) {
2917 m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][0];
2918 m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][2];
2919 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
6be22b3f 2920 }
2921 }
2922 //
2923 m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
2924 } while( (fCurrentModule=fCurrentModule->GetParent()) );
2925 //
2926 if (nmodTested>0 && zeroX) {AliInfo("Skipping: zero global X derivatives!");return -1;}
2927 if (nmodTested>0 && zeroZ) {AliInfo("Skipping: zero global Z derivatives!");return -1;}
2928 //
2929 // ok, can copy to m
2930 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2931 m.fMeas[kX] = fMeasLoc[0]-fPintLoc0[0];
2932 m.fSigma[kX] = fSigmaLoc[0];
2933 //
2934 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2935 m.fMeas[kZ] = fMeasLoc[2]-fPintLoc0[2];
2936 m.fSigma[kZ] = fSigmaLoc[2];
2937 //
2938 m.fNGlobFilled = ifill;
2939 fCurrentModule = endModule;
2940 //
8fd71c0a 2941 status += Int_t(!zeroX && !zeroZ); // 0 - only locals, 1 locals + globals
2942 return status;
6be22b3f 2943}
2944
2945//________________________________________________________________________________________________________
2946Int_t AliITSAlignMille2::AddLocalEquationTPA(Mille2Data &m)
2947{
2948 /// Define local equation for current cluster in X Y and Z coor.
2949 /// and store them to memory
2950 /// return -1 in case of failure to build some equation
2951 /// 0 if no free global parameters were found but local eq is built
2952 /// 1 if both local and global eqs are built
2953 //
8102b2c9 2954 static int cnt = 0;
2955 Bool_t dbg = kFALSE;//kTRUE;
2956 if (++cnt>100000) dbg = kFALSE;
2957
6be22b3f 2958 int curpoint = fCluster.GetUniqueID()-1;
2959 TGeoHMatrix *tempHMat = GetSensorCurrMatrixSID(fCurrentSensID);// fCurrentModule->GetSensitiveVolumeMatrix(fCluster.GetVolumeID());
2960 //
2961 fTPAFitter->GetDResDParams(&fDerivativeLoc[0][0], curpoint); // resid. derivatives over the track parameters
8693b6b7 2962 if (fCurvFitWasConstrained && fFixCurvIfConstraned && !IsZero(fBField))
2963 for (int i=3;i--;) fDerivativeLoc[AliITSTPArrayFit::kR0][i] = 0; //Fix curvarute
2964 //
6be22b3f 2965 for (Int_t i=fNLocal; i--;) tempHMat->MasterToLocalVect(fDerivativeLoc[i],m.fDerLoc[i]);
2966 //
8fd71c0a 2967 int status = 0;
6be22b3f 2968 // derivatives over the global parameters ---------------------------------------->>>
ef24eb3b 2969 Double_t dGL[3]; // derivative of global position vs local X (for SDD)
6be22b3f 2970 Double_t dRdP[3][3]; // derivative of local residuals vs local position
2971 Double_t dPdG[AliITSAlignMille2Module::kMaxParGeom][3]; // derivatives of local position vs global params
2972 fTPAFitter->GetDResDPos(&fDerivativeGlo[0][0], curpoint);
8102b2c9 2973 if (fCurrentSensID!=kVtxSensID) for (int i=3;i--;) tempHMat->MasterToLocalVect(fDerivativeGlo[i],dRdP[i]);
2974 else for (int i=3;i--;) for (int j=3;j--;) dRdP[i][j] = fDerivativeGlo[i][j]; // vertex constraint is in lab
2975 //
2976 if (dbg) {
2977 printf("\nCurrentMod: %s Sens:%d\n",fCurrentModule->GetName(),fCurrentSensID); //RRR
2978 printf("Module Matrix: ");
2979 fCurrentModule->GetMatrix()->Print(); //RRR
2980 for (int i=0;i<3;i++) {
2981 printf("dRdP[M%d][resI] ",i); for (int j=0;j<3;j++) printf(":[%d] %+.3e ",j,dRdP[i][j]); printf("\n");
2982 }//RRR
2983 printf("Sensor Matrix: "); tempHMat->Print();
2984 }
6be22b3f 2985 UInt_t ifill=0, dfDone = 0;
2986 m.fNModFilled = 0;
2987 //
2988 AliITSAlignMille2Module* endModule = fCurrentModule;
2989 //
8102b2c9 2990 m.fModuleID[0] = fCurrentModule->GetUniqueID(); // always register id of the base module, even if it has no DOF
2991 //
6be22b3f 2992 do {
2993 if (fCurrentModule->GetNParFree()==0) continue;
8fd71c0a 2994 status = 1;
6be22b3f 2995 if (!fUseGlobalDelta) dfDone = 0; // for local deltas the derivatives at diff. levels are different
2996 Bool_t jacobOK = kFALSE;
2997 //
2998 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2999 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
3000 //
3001 if (!TestWordBit(dfDone,i)) { // need to calculate new derivative
ef24eb3b 3002 if (!jacobOK) {
8102b2c9 3003 if (fCurrentSensID!=kVtxSensID) {
3004 fCurrentModule->CalcDerivDPosDPar(fCluster.GetVolumeID(),fMeasLoc,&dPdG[0][0]);
3005 if (dbg) {
3006 for (int i1=0;i1<3;i1++) {
3007 printf("Jacob:dPdG[gpar%d][Mj]",i1); for (int j1=0;j1<3;j1++) printf(":[%d] %+.3e ",j1,dPdG[i1][j1]); printf("\n");//RRR
3008 }
3009 }
3010 }
3011 else {
3012 // this is a vertex constraint: only lateral shifts are allowed, no rotations
3013 for (int ip=AliITSAlignMille2Module::kMaxParGeom;ip--;) for (int jp=3;jp--;) dPdG[ip][jp] = (ip==jp) ? 1:0;
3014 }
ef24eb3b 3015 jacobOK = kTRUE;
3016 }
6be22b3f 3017 // dRes_j/dGlo_i = \sum_{k=1:3} dRes_j/dPos_k * dPos_k/dGlo_i
3018 fDerivativeGlo[i][kX] = dRdP[kX][kX]*dPdG[i][kX] + dRdP[kY][kX]*dPdG[i][kY] + dRdP[kZ][kX]*dPdG[i][kZ];
3019 fDerivativeGlo[i][kY] = dRdP[kX][kY]*dPdG[i][kX] + dRdP[kY][kY]*dPdG[i][kY] + dRdP[kZ][kY]*dPdG[i][kZ];
3020 fDerivativeGlo[i][kZ] = dRdP[kX][kZ]*dPdG[i][kX] + dRdP[kY][kZ]*dPdG[i][kY] + dRdP[kZ][kZ]*dPdG[i][kZ];
3021 SetWordBit(dfDone,i);
3022 }
3023 //
8102b2c9 3024 if (dbg) {
3025 printf("Level %s DGlob[par%d][resJ] ",fCurrentModule->GetName(),i); //RRR
3026 for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k, fDerivativeGlo[i][k]); printf("\n");//RRR
3027 }
6be22b3f 3028 m.fDerGlo[ifill][kX] = fDerivativeGlo[i][kX];
3029 m.fDerGlo[ifill][kY] = fDerivativeGlo[i][kY];
3030 m.fDerGlo[ifill][kZ] = fDerivativeGlo[i][kZ];
3031 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(i);
3032 //
3033 }
3034 //
3035 if ( fCurrentModule->IsSDD() ) { // specific for SDD
3036 //
3037 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
3038 // where V0 and T are the nominal drift velocity, time and time0
3039 // and the dT0 and dV are the corrections:
ef24eb3b 3040 // drloc_i/dT0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dT0 =
3041 // = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dT0
3042 // = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * V0
3043 //
3044 // drloc_i/dV0 = sum_j drloc_i/dMeasGlo_j * dMeasGlo_j/dV0 =
3045 // = sum_j drloc_i/dMeasGlo_j sum_k dMeasGlo_j/dMeasLoc_k * dMeasLoc_k/dV0
3046 // = sum_j drloc_i/dMeasGlo_j dMeasGlo_j/dMeasLoc_X * T0
3047
6be22b3f 3048 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
3049 //
ef24eb3b 3050 Bool_t jacOK = kFALSE;
3051 //Int_t sddLR = fMeasLoc[kX]>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
3052 Int_t sddLR = GetVDriftSDD()>0 ? AliITSAlignMille2Module::kDOFDVL : AliITSAlignMille2Module::kDOFDVR;
6be22b3f 3053 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
3054 if (!TestWordBit(dfDone, AliITSAlignMille2Module::kDOFT0)) {
3055 double vdrift = GetVDriftSDD();
ef24eb3b 3056 JacobianPosGloLoc(kX,dGL);
3057 jacOK = kTRUE;
3058 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX] =
3059 vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3060 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY] =
3061 vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3062 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ] =
3063 vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3064 //
6be22b3f 3065 SetWordBit(dfDone, AliITSAlignMille2Module::kDOFT0);
3066 }
3067 m.fDerGlo[ifill][kX] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kX];
3068 m.fDerGlo[ifill][kY] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kY];
3069 m.fDerGlo[ifill][kZ] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][kZ];
3070 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
3071 }
3072 //
ef24eb3b 3073 if (fCurrentModule->GetParOffset(sddLR)>=0) {
3074 if (!TestWordBit(dfDone, sddLR)) {
6be22b3f 3075 double tdrift = TMath::Sign(GetTDriftSDD(), GetVDriftSDD());
ef24eb3b 3076 double vdrift = fIsSDDVDriftMult ? TMath::Abs(GetVDriftSDD()) : 1;
3077 if (!jacOK) JacobianPosGloLoc(kX,dGL);
3078 fDerivativeGlo[sddLR][kX] =
3079 -tdrift*vdrift*(dRdP[kX][kX]*dGL[kX] + dRdP[kY][kX]*dGL[kY] + dRdP[kZ][kX]*dGL[kZ]);
3080 fDerivativeGlo[sddLR][kY] =
3081 -tdrift*vdrift*(dRdP[kX][kY]*dGL[kX] + dRdP[kY][kY]*dGL[kY] + dRdP[kZ][kY]*dGL[kZ]);
3082 fDerivativeGlo[sddLR][kZ] =
3083 -tdrift*vdrift*(dRdP[kX][kZ]*dGL[kX] + dRdP[kY][kZ]*dGL[kY] + dRdP[kZ][kZ]*dGL[kZ]);
3084 SetWordBit(dfDone, sddLR);
6be22b3f 3085 }
ef24eb3b 3086 m.fDerGlo[ifill][kX] = fDerivativeGlo[sddLR][kX];
3087 m.fDerGlo[ifill][kY] = fDerivativeGlo[sddLR][kY];
3088 m.fDerGlo[ifill][kZ] = fDerivativeGlo[sddLR][kZ];
3089 m.fParMilleID[ifill++] = fCurrentModule->GetParOffset(sddLR);
6be22b3f 3090 }
3091 }
3092 //
3093 m.fModuleID[m.fNModFilled++] = fCurrentModule->GetUniqueID();
3094 } while( (fCurrentModule=fCurrentModule->GetParent()) );
3095 //
3096 // store first local residuals
3097 fTPAFitter->GetResiduals(fPintLoc , curpoint); // lab residuals
3098 for (int i=3;i--;) fPintLoc[i] = -fPintLoc[i];
8102b2c9 3099 if (fCurrentSensID!=kVtxSensID) tempHMat->MasterToLocalVect(fPintLoc,m.fMeas); // local residuals
3100 else for (int i=3;i--;) m.fMeas[i] = fPintLoc[i];
3101 if (dbg) {
3102 printf("res(meas-loc) "); for (int k=0;k<3;k++) printf(":[%d] %+.3e ",k,m.fMeas[k]); printf("\n");
3103 printf("Fin:%s %+e %+e\n",endModule->GetName(), fDerivativeGlo[kZ][kZ], fPintLoc[kZ]);
3104 }//RRR
6be22b3f 3105 m.fSigma[kX] = fSigmaLoc[kX];
3106 m.fSigma[kY] = fSigmaLoc[kY];
3107 m.fSigma[kZ] = fSigmaLoc[kZ];
3108 //
3109 m.fNGlobFilled = ifill;
3110 fCurrentModule = endModule;
3111 //
8fd71c0a 3112 return status;
6be22b3f 3113}
3114
3115//________________________________________________________________________________________________________
3116void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
3117{
3118 /// Set local equations with data stored in m
3119 /// return 0 if success
3120 //
8693b6b7 3121 Bool_t locPatt[kNLocal] = {0}; // pattern of lacal eq's to account
3122 for (int i=fNLocal; i--;) {
3123 if (locPatt[i]) continue; // already set
3124 for (Int_t j=0; j<neq; j++) for (int ic=3;ic--;) if (!IsZero(marr[j].fDerLoc[i][ic])) locPatt[i]=kTRUE;
3125 }
3126 //
6be22b3f 3127 for (Int_t j=0; j<neq; j++) {
3128 //
3129 const Mille2Data &m = marr[j];
3130 //
3131 Bool_t filled = kFALSE;
3132 for (int ic=3;ic--;) {
ef24eb3b 3133 // for the diamond point (if any) the Y residual is accounted
3134 if (ic==kY && !fUseLocalYErr && !(m.fModuleID[0]==fDiamondModID)) continue;
6be22b3f 3135 AliDebug(2,Form("setting local equation %c with fMeas=%.6f and fSigma=%.6f",fgkXYZ[ic],m.fMeas[ic], m.fSigma[ic]));
8693b6b7 3136 Int_t nzero = 0, naddl = 0;
3137 for (int i=0;i<=fNLocal;i++) if (locPatt[i]) nzero += SetLocalDerivative(naddl++,m.fDerLoc[i][ic] );
ef24eb3b 3138 if (nzero==fNLocal) {
3139 AliInfo(Form("Skipping %c residual due to the zero derivatives!",fgkXYZ[ic]));
3140 continue;
3141 }
8fd71c0a 3142 for (int i=m.fNGlobFilled;i--;) SetGlobalDerivative( m.fParMilleID[i] , m.fDerGlo[i][ic] );
6be22b3f 3143 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.fMeas[ic], m.fSigma[ic]);
3144 filled = kTRUE;
3145 //
3146 }
3147 //
3148 if (filled) for (int i=m.fNModFilled;i--;) GetMilleModule(m.fModuleID[i])->IncNProcessedPoints();
3149 }
ef24eb3b 3150 //
3151 double wgh = 1;
3152 if (GetWeightPt() && fTPAFitter) {
3153 wgh = fTPAFitter->GetPt();
3154 if (wgh>10) wgh = 10.;
3155 if (wgh<0) wgh = fTPAFitter->IsTypeCosmics() ? 7 : 0.5;
3156 if (GetWeightPt()>0) wgh = TMath::Power(wgh,GetWeightPt());
3157 }
3158 fMillepede->SetRecordWeight(wgh*fTrackWeight);
8102b2c9 3159 fMillepede->SetRecordRun(fRunID);
ef24eb3b 3160 //
6be22b3f 3161}
3162
3163//________________________________________________________________________________________________________
3164Int_t AliITSAlignMille2::GlobalFit()
3165{
3166 /// Call global fit; Global parameters are stored in parameters
3167 if (!fIsMilleInit) Init();
3168 //
3169 ApplyPreConstraints();
3170 int res = fMillepede->GlobalFit();
3171 AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
3172 if (res) {
3173 // fetch the parameters
3174 for (int imd=fNModules;imd--;) {
3175 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3176 int nprocp = 0;
3177 for (int ip=mod->GetNParTot();ip--;) {
3178 int idp = mod->GetParOffset(ip);
3179 if (idp<0) continue; // was not in the explicit fit
3180 mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
3181 mod->SetParErr(ip,fMillepede->GetFinalError(idp));
3182 int np = fMillepede->GetProcessedPoints(idp);
3183 if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
3184 }
3185 if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
3186 }
3187
3188 }
3189 ApplyPostConstraints();
3190 return res;
3191}
3192
3193//________________________________________________________________________________________________________
3194void AliITSAlignMille2::PrintGlobalParameters()
3195{
3196 /// Print global parameters
3197 if (!fIsMilleInit) {
3198 AliInfo("Millepede has not been initialized!");
3199 return;
3200 }
3201 fMillepede->PrintGlobalParameters();
3202}
3203
3204//________________________________________________________________________________________________________
3205Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
3206{
3207 // load definitions of supermodules from a root file
3208 // return 0 if success
ef24eb3b 3209 AliInfo(Form("Loading SuperModule definitions from %s",sfile));
6be22b3f 3210 TFile *smf=TFile::Open(sfile);
3211 if (!smf->IsOpen()) {
3212 AliInfo(Form("Cannot open supermodule file %s",sfile));
3213 return -1;
3214 }
3215
3216 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
3217 if (!sma) {
3218 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
3219 return -2;
3220 }
3221 Int_t nsma=sma->GetEntriesFast();
3222 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
3223 //
3224 // pepo200709
3225 Char_t st[2048];
3226 char symname[250];
3227 // end pepo200709
3228
3229 UShort_t volid;
3230 TGeoHMatrix m;
3231 //
3232 for (Int_t i=0; i<nsma; i++) {
3233 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
3234 volid=a->GetVolUID();
c82cdb65 3235 strncpy(st,a->GetSymName(),TMath::Min(sizeof(st),strlen(a->GetSymName())+1));
6be22b3f 3236 a->GetMatrix(m);
3237 //
b80c197e 3238 memset(symname,0,250*sizeof(char));
c82cdb65 3239 sscanf(st,"%249s",symname);
6be22b3f 3240 //
3241 // decode module list
3242 char *stp=strstr(st,"ModuleList:");
3243 if (!stp) return -3;
3244 stp += 11;
3245 int idx[2200];
3246 char spp[200]; int jp=0;
3247 char cl[20];
c82cdb65 3248 strncpy(st,stp,TMath::Min(sizeof(st),strlen(stp)+1));
6be22b3f 3249 int l=strlen(st);
3250 int j=0;
3251 int n=0;
3252 //
3253 while (j<=l) {
3254 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
3255 spp[jp]=0;
3256 jp=0;
3257 if (strlen(spp)) {
3258 int k=strcspn(spp,"-");
3259 if (k<int(strlen(spp))) { // c'e' il -
c82cdb65 3260 strncpy(cl,&(spp[k+1]), TMath::Min(sizeof(cl),strlen(&spp[k+1])+1));
6be22b3f 3261 spp[k]=0;
3262 int ifrom=atoi(spp); int ito=atoi(cl);
3263 for (int b=ifrom; b<=ito; b++) {
3264 idx[n]=b;
3265 n++;
3266 }
3267 }
3268 else { // numerillo singolo
3269 idx[n]=atoi(spp);
3270 n++;
3271 }
3272 }
3273 }
3274 else {
3275 spp[jp]=st[j];
3276 jp++;
3277 }
3278 j++;
3279 }
3280 UShort_t volidsv[2198];
3281 for (j=0;j<n;j++) {
3282 volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
3283 if (!volidsv[j]) {
3284 AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
3285 return -5;
3286 }
3287 }
3288 Int_t smindex=int(2198+volid-14336); // virtual index
3289 //
3290 fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
3291 //
3292 fNSuperModules++;
3293 }
3294 //
3295 smf->Close();
3296 //
3297 return 0;
3298}
3299
3300//________________________________________________________________________________________________________
3301void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
3302{
3303 // require that sum of modifications for the childs of this module is = val, i.e.
3304 // the internal corrections moves the module as a whole by fixed value (0 by default).
3305 // pattern is the bit pattern for the parameters to constrain
3306 //
3307 if (fIsMilleInit) {
3308 AliInfo("Millepede has been already initialized: no constrain may be added!");
3309 return;
3310 }
3311 if (!GetMilleModule(idm)->GetNChildren()) return;
3312 TString nm = "cstrSUMean";
3313 nm += GetNConstraints();
3314 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3315 idm,val,pattern);
3316 cstr->SetConstraintID(GetNConstraints());
3317 fConstraints.Add(cstr);
3318}
3319
3320//________________________________________________________________________________________________________
3321void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
3322{
3323 // require that median of the modifications for the childs of this module is = val, i.e.
3324 // the internal corrections moves the module as a whole by fixed value (0 by default)
3325 // module the outliers.
3326 // pattern is the bit pattern for the parameters to constrain
3327 // The difference between the mean and the median will be transfered to the parent
3328 if (fIsMilleInit) {
3329 AliInfo("Millepede has been already initialized: no constrain may be added!");
3330 return;
3331 }
3332 if (!GetMilleModule(idm)->GetNChildren()) return;
3333 TString nm = "cstrSUMed";
3334 nm += GetNConstraints();
3335 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3336 idm,val,pattern);
3337 cstr->SetConstraintID(GetNConstraints());
3338 fConstraints.Add(cstr);
3339}
3340
3341//________________________________________________________________________________________________________
3342void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
3343{
3344 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3345 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3346 // pattern is the bit pattern for the parameters to constrain
3347 //
3348 if (fIsMilleInit) {
3349 AliInfo("Millepede has been already initialized: no constrain may be added!");
3350 return;
3351 }
3352 TString nm = "cstrOMean";
3353 nm += GetNConstraints();
3354 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
3355 -1,val,pattern);
3356 cstr->SetConstraintID(GetNConstraints());
3357 fConstraints.Add(cstr);
3358}
3359
3360//________________________________________________________________________________________________________
3361void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
3362{
3363 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3364 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3365 // pattern is the bit pattern for the parameters to constrain
3366 //
3367 if (fIsMilleInit) {
3368 AliInfo("Millepede has been already initialized: no constrain may be added!");
3369 return;
3370 }
3371 TString nm = "cstrOMed";
3372 nm += GetNConstraints();
3373 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
3374 -1,val,pattern);
3375 cstr->SetConstraintID(GetNConstraints());
3376 fConstraints.Add(cstr);
3377}
3378
3379//________________________________________________________________________________________________________
3380void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
3381{
3382 // apply constraint on parameters in the local frame
3383 if (fIsMilleInit) {
3384 AliInfo("Millepede has been already initialized: no constrain may be added!");
3385 return;
3386 }
3387 AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
3388 cstr->SetConstraintID(GetNConstraints());
3389 fConstraints.Add(cstr);
3390}
3391
3392//________________________________________________________________________________________________________
3393void AliITSAlignMille2::ApplyGaussianConstraint(const AliITSAlignMille2ConstrArray* cstr)
3394{
3395 // apply the constraint on the local corrections of a list of modules
3396 int nmod = cstr->GetNModules();
3397 double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
3398 //
ef24eb3b 3399 // check if this not special SDDT0 constraint
3400 if (cstr->GetPattern()==BIT(AliITSAlignMille2Module::kDOFT0)) {
3401 for (int i=0;i<cstr->GetNModules()-1;i++) {
3402 AliITSAlignMille2Module *mdI = GetMilleModule(cstr->GetModuleID(i));
3403 if (!mdI->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3404 for (int j=i+1;j<cstr->GetNModules();j++) {
3405 AliITSAlignMille2Module *mdJ = GetMilleModule(cstr->GetModuleID(j));
3406 if (!mdJ->IsFreeDOF(AliITSAlignMille2Module::kDOFT0)) continue;
3407 //
3408 ResetLocalEquation();
3409 fGlobalDerivatives[mdI->GetParOffset(AliITSAlignMille2Module::kDOFT0)] = 1;
3410 fGlobalDerivatives[mdJ->GetParOffset(AliITSAlignMille2Module::kDOFT0)] =-1;
3411 AddConstraint(fGlobalDerivatives, 0, 1.E-6);
3412 }
3413 }
3414 return;
3415 }
3416
6be22b3f 3417 for (int imd=nmod;imd--;) {
3418 int modID = cstr->GetModuleID(imd);
3419 AliITSAlignMille2Module* mod = GetMilleModule(modID);
3420 ResetLocalEquation();
3421 int nadded = 0;
3422 double value = cstr->GetValue();
3423 double sigma = cstr->GetError();
3424 //
3425 // in case the reference (survey) deltas were imposed for Gaussian constraints
3426 // already accumulated corrections: they must be subtracted from the constraint value.
3427 if (IsConstraintWrtRef()) {
3428 //
3429 Double_t precal[AliITSAlignMille2Module::kMaxParTot];
3430 Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
3431 for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
3432 //
3433 // check if there was a reference delta provided for this module
3434 AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
3435 if (parref) parref->GetPars(refcal, refcal+3); // found reference delta
3436 //
3437 // extract already applied local corrections for this module
3438 if (fPrealignment) {
3439 //
3440 AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
3441 if (preo) {
3442 TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); // Delta_Glob * Delta_Glob_Par * M
3443 preo->GetMatrix(preMat); // Delta_Glob
3444 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
3445 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
3446 AliAlignObjParams algob;
3447 algob.SetMatrix(tmpMat);
3448 algob.GetPars(precal,precal+3); // local corrections for geometry
3449 }
3450 }
3451 //
3452 // subtract the contribution to constraint from precalibration
3453 for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
3454 //
3455 }
3456 //
3457 if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
3458 //
3459 for (int ipar=cstr->GetNCoeffs();ipar--;) {
3460 double coef = cstr->GetCoeff(ipar);
3461 if (IsZero(coef)) continue;
3462 //
3463 if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
3464 // we are working with local params or if the given param is not related to geometry,
3465 // apply the constraint directly
3466 int parPos = mod->GetParOffset(ipar);
3467 if (parPos<0) continue; // not in the fit
3468 fGlobalDerivatives[parPos] += coef;
3469 nadded++;
3470 }
3471 else { // we are working with global params, while the constraint is on local ones -> jacobian
3472 for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
3473 int parPos = mod->GetParOffset(jpar);
3474 if (parPos<0) continue;
3475 fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
3476 nadded++;
3477 }
3478 }
3479 }
3480 if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
3481 }
3482 //
3483}
3484
3485//________________________________________________________________________________________________________
3486void AliITSAlignMille2::ApplyPreConstraints()
3487{
3488 // apply constriants which cannot be imposed after the fit
3489 int nconstr = GetNConstraints();
3490 for (int i=0;i<nconstr;i++) {
3491 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3492 //
3493 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
3494 ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
3495 continue;
3496 }
3497 //
3498 if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
3499 //
3500 if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
3501 // apply constraint on the mean's before the fit
3502 int imd = cstr->GetModuleID();
3503 if (imd>=0) {
3504 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3505 UInt_t pattern = 0;
3506 for (int ipar=mod->GetNParTot();ipar--;) {
3507 if (!cstr->IncludesParam(ipar)) continue;
3508 if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
3509 pattern |= 0x1<<ipar;
3510 cstr->SetApplied(ipar);
3511 }
3512 ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
3513 //
3514 }
3515 else if (!PseudoParentsAllowed()) {
3516 ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
3517 cstr->SetApplied(-1);
3518 }
3519 }
ef24eb3b 3520 //
3521 // do we need to tie the SDD left/right VDrift corrections
3522 for (int imd=0;imd<fNModules;imd++) {
3523 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3524 if (mod->IsSDD() && mod->IsVDriftLRSame()) TieSDDVDriftsLR(mod);
3525 }
3526 //
6be22b3f 3527}
3528
3529//________________________________________________________________________________________________________
3530void AliITSAlignMille2::ApplyPostConstraints()
3531{
3532 // apply constraints which can be imposed after the fit
3533 int nconstr = GetNConstraints();
3534 Bool_t convGlo = kFALSE;
3535 // check if there is something to do
3536 int ntodo = 0;
3537 for (int i=0;i<nconstr;i++) {
3538 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3539 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3540 if (cstr->GetRemainingPattern() == 0) continue;
3541 ntodo++;
3542 }
3543 if (!ntodo) return;
3544 //
3545 if (!fUseGlobalDelta) { // need to convert to global params
3546 ConvertParamsToGlobal();
3547 convGlo = kTRUE;
3548 }
3549 //
3550 for (int i=0;i<nconstr;i++) {
3551 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
3552 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
3553 //
3554 int imd = cstr->GetModuleID();
3555 //
3556 if (imd>=0) {
3557 AliITSAlignMille2Module* mod = GetMilleModule(imd);
8102b2c9 3558 if (mod->IsNotInConf()) continue;
6be22b3f 3559 UInt_t pattern = 0;
3560 for (int ipar=mod->GetNParTot();ipar--;) {
3561 if (cstr->IsApplied(ipar)) continue;
3562 if (!cstr->IncludesParam(ipar)) continue;
3563 if (!mod->IsFreeDOF(ipar)) continue; // parameter is fixed, will not apply constraint
3564 pattern |= 0x1<<ipar;
3565 cstr->SetApplied(ipar);
3566 }
3567 if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
3568 //
3569 }
3570 else if (PseudoParentsAllowed()) {
3571 UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
3572 PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
3573 cstr->SetApplied(-1);
3574 }
3575 }
3576 // if there was a conversion, rewind it
3577 if (convGlo) ConvertParamsToLocal();
3578 //
3579}
3580
3581//________________________________________________________________________________________________________
3582void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
3583{
3584 // require that sum of modifications for the childs of this module is = val, i.e.
3585 // the internal corrections moves the module as a whole by fixed value (0 by default).
3586 // pattern is the bit pattern for the parameters to constrain
3587 //
3588 //
3589 AliITSAlignMille2Module* mod = GetMilleModule(idm);
3590 //
3591 for (int ip=0;ip<kNParCh;ip++) {
3592 if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
3593 ResetLocalEquation();
3594 int nadd = 0;
3595 for (int ich=mod->GetNChildren();ich--;) {
3596 int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
3597 if (idpar<0) continue;
3598 fGlobalDerivatives[idpar] = 1.0;
3599 nadd++;
3600 }
3601 //
3602 if (nadd>0) {
3603 AddConstraint(fGlobalDerivatives,val);
3604 AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
3605 }
3606 }
3607 //
3608}
3609
3610//________________________________________________________________________________________________________
3611void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
3612{
3613 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
3614 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
3615 // pattern is the bit pattern for the parameters to constrain
3616 //
3617 for (int ip=0;ip<kNParCh;ip++) {
3618 //
3619 if ( !((pattern>>ip)&0x1) ) continue;
3620 ResetLocalEquation();
3621 int nadd = 0;
3622 for (int imd=fNModules;imd--;) {
3623 AliITSAlignMille2Module* mod = GetMilleModule(imd);
8102b2c9 3624 if (mod->IsNotInConf()) continue; // dummy module
3625 AliITSAlignMille2Module* par = mod->GetParent();
3626 while (par && par->IsNotInConf() ) par = par->GetParent(); // use only decalred parents
3627 if (par) continue; // this is not an orphan
6be22b3f 3628 int idpar = mod->GetParOffset(ip);
3629 if (idpar<0) continue;
3630 fGlobalDerivatives[idpar] = 1.0;
3631 nadd++;
3632 }
3633 if (nadd>0) {
3634 AddConstraint(fGlobalDerivatives,val);
3635 AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
3636 }
3637 }
3638 //
3639 //
3640}
3641
3642//________________________________________________________________________________________________________
3643void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
3644{
3645 // require that median or mean of the modifications for the childs of this module is = val, i.e.
3646 // the internal corrections moves the module as a whole by fixed value (0 by default)
3647 // module the outliers.
3648 // pattern is the bit pattern for the parameters to constrain
3649 // The difference between the mean and the median will be transfered to the parent
3650 //
3651 AliITSAlignMille2Module* parent = GetMilleModule(idm);
3652 int nc = parent->GetNChildren();
3653 //
3654 double *tmpArr = new double[nc];
3655 //
3656 for (int ip=0;ip<kNParCh;ip++) {
3657 int npc = 0;
3658 if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
3659 // compute the mean and median of the deltas
3660 int nfree = 0;
3661 for (int ich=nc;ich--;) {
3662 AliITSAlignMille2Module* child = parent->GetChild(ich);
3663 // if (!child->IsFreeDOF(ip)) continue;
3664 tmpArr[nfree++] = child->GetParVal(ip);
3665 }
3666 double median=0,mean=0;
3667 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3668 mean += tmpArr[ic0];
3669 for (int ic1=ic0+1;ic1<nfree;ic1++)
3670 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3671 }
3672 //
3673 int kmed = nfree/2;
3674 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3675 if (nfree>0) mean /= nfree;
3676 //
3677 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3678 //
3679 for (int ich=nc;ich--;) {
3680 AliITSAlignMille2Module* child = parent->GetChild(ich);
3681 // if (!child->IsFreeDOF(ip)) continue;
3682 child->SetParVal(ip, child->GetParVal(ip) + shift);
3683 npc++;
3684 }
3685 //
3686 parent->SetParVal(ip, parent->GetParVal(ip) - shift);
ef24eb3b 3687 AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d children of module %d: %s",
6be22b3f 3688 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3689 ip,npc,idm,parent->GetName()));
3690 }
3691 delete[] tmpArr;
3692 //
3693 //
3694}
3695
3696//________________________________________________________________________________________________________
3697void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
3698{
3699 // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
3700 // the corrections moves the whole setup by fixed value (0 by default).
3701 // pattern is the bit pattern for the parameters to constrain
3702 //
3703 int nc = fNModules;
3704 //
3705 int norph = 0;
8102b2c9 3706 for (int ich=nc;ich--;) {
3707 AliITSAlignMille2Module *par= GetMilleModule(ich)->GetParent();
3708 while (par && par->IsNotInConf()) par = par->GetParent(); // use only decalred parents
3709 if (!par) norph ++;
3710 }
3711 //
6be22b3f 3712 if (!norph) return;
3713 double *tmpArr = new double[norph];
8102b2c9 3714 for (int i=norph;i--;) tmpArr[i] = 0;
6be22b3f 3715 //
3716 for (int ip=0;ip<kNParCh;ip++) {
3717 int npc = 0;
3718 if ( !((pattern>>ip)&0x1)) continue;
3719 // compute the mean and median of the deltas
3720 int nfree = 0;
3721 for (int ich=nc;ich--;) {
3722 AliITSAlignMille2Module* child = GetMilleModule(ich);
8102b2c9 3723 if (child->IsNotInConf()) continue; // dummy module
6be22b3f 3724 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
8102b2c9 3725 AliITSAlignMille2Module* par = child->GetParent();
3726 while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3727 if (par) continue;
6be22b3f 3728 tmpArr[nfree++] = child->GetParVal(ip);
3729 }
3730 double median=0,mean=0;
3731 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
3732 mean += tmpArr[ic0];
3733 for (int ic1=ic0+1;ic1<nfree;ic1++)
3734 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
3735 }
3736 //
3737 int kmed = nfree/2;
3738 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
3739 if (nfree>0) mean /= nfree;
3740 //
3741 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
3742 //
3743 for (int ich=nc;ich--;) {
3744 AliITSAlignMille2Module* child = GetMilleModule(ich);
8102b2c9 3745 if (child->IsNotInConf()) continue; // dummy module
6be22b3f 3746 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
8102b2c9 3747 AliITSAlignMille2Module* par = child->GetParent();
3748 while (par && par->IsNotInConf()) par = par->GetParent(); // count only declared parents
3749 if (par) continue;
6be22b3f 3750 child->SetParVal(ip, child->GetParVal(ip) + shift);
3751 npc++;
3752 }
3753 //
ef24eb3b 3754 AliInfo(Form("%s constraint: added %+f shift to param[%d] of %d orphan modules",
6be22b3f 3755 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
3756 ip,npc));
3757 }
3758 delete[] tmpArr;
3759 //
3760}
3761
3762//________________________________________________________________________________________________________
3763Bool_t AliITSAlignMille2::IsParModConstrained(const AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
3764{
3765 // check if par of the module participates in some constraint, and set the flag for their types
3766 meanmed = gaussian = kFALSE;
3767 //
3768 if ( mod->IsParConstrained(par) ) gaussian = kTRUE; // direct constraint on this param
3769 //
3770 for (int icstr=GetNConstraints();icstr--;) {
3771 AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
3772 //
3773 if (!cstr->IncludesModPar(mod,par)) continue;
3774 if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
3775 else meanmed = kTRUE;
3776 //
3777 if (meanmed && gaussian) break; // no sense to check further
3778 }
3779 //
3780 return meanmed||gaussian;
3781}
3782
3783//________________________________________________________________________________________________________
3784Bool_t AliITSAlignMille2::IsParModFamilyVaried(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3785{
3786 // check if parameter par is varied for this module or its children up to the level depth
3787 if (depth<0) return kFALSE;
3788 if (mod->GetParOffset(par)>=0) return kTRUE;
3789 for (int icld=mod->GetNChildren();icld--;) {
3790 AliITSAlignMille2Module* child = mod->GetChild(icld);
3791 if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
3792 }
3793 return kFALSE;
3794 //
3795}
3796
3797/*
3798//________________________________________________________________________________________________________
3799Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3800{
3801 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3802 if (depth<0) return kTRUE;
3803 for (int icld=mod->GetNChildren();icld--;) {
3804 AliITSAlignMille2Module* child = mod->GetChild(icld);
3805 //if (child->GetParOffset(par)<0) continue; // fixed
3806 Bool_t cstMM=kFALSE,cstGS=kFALSE;
3807 // does this child have gaussian constraint ?
3808 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3809 // check its children
3810 if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
3811 }
3812 return kFALSE;
3813 //
3814}
3815*/
3816
3817//________________________________________________________________________________________________________
3818Bool_t AliITSAlignMille2::IsParFamilyFree(const AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
3819{
3820 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
3821 if (depth<0) return kFALSE;
3822 for (int icld=mod->GetNChildren();icld--;) {
3823 AliITSAlignMille2Module* child = mod->GetChild(icld);
3824 //if (child->GetParOffset(par)<0) continue; // fixed
3825 Bool_t cstMM=kFALSE,cstGS=kFALSE;
3826 // does this child have gaussian constraint ?
3827 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
3828 // check its children
3829 if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
3830 }
3831 return kFALSE;
3832 //
3833}
3834
3835//________________________________________________________________________________________________________
3836Double_t AliITSAlignMille2::GetTDriftSDD() const
3837{
3838 // obtain drift time corrected for t0
3839 double t = fCluster.GetDriftTime();
3840 return t - fDriftTime0[ fCluster.GetUniqueID()-1 ];
3841}
3842
3843//________________________________________________________________________________________________________
3844Double_t AliITSAlignMille2::GetVDriftSDD() const
3845{
3846 // obtain corrected drift speed
3847 return fDriftSpeed[ fCluster.GetUniqueID()-1 ];
3848}
3849
3850//________________________________________________________________________________________________________
3851Bool_t AliITSAlignMille2::FixedOrphans() const
3852{
3853 // are there fixed modules with no parent (normally in such a case
3854 // the constraints on the orphans should not be applied
3855 if (!IsConfigured()) {
3856 AliInfo("Still not configured");
3857 return kFALSE;
3858 }
3859 for (int i=0;i<fNModules;i++) {
3860 AliITSAlignMille2Module* md = GetMilleModule(i);
8102b2c9 3861 if (md->IsNotInConf()) continue;
6be22b3f 3862 if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
3863 }
3864 return kFALSE;
3865}
3866
3867//________________________________________________________________________________________________________
a878a855 3868void AliITSAlignMille2::ConvertParamsToGlobal() const
6be22b3f 3869{
3870 // convert params in local frame to global one
3871 double pars[AliITSAlignMille2Module::kMaxParGeom];
3872 for (int imd=fNModules;imd--;) {
3873 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3874 if (mod->GeomParamsGlobal()) continue;
3875 mod->GetGeomParamsGlo(pars);
3876 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3877 mod->SetGeomParamsGlobal(kTRUE);
3878 }
3879}
3880
3881//________________________________________________________________________________________________________
a878a855 3882void AliITSAlignMille2::ConvertParamsToLocal() const
6be22b3f 3883{
3884 // convert params in global frame to local one
3885 double pars[AliITSAlignMille2Module::kMaxParGeom];
3886 for (int imd=fNModules;imd--;) {
3887 AliITSAlignMille2Module* mod = GetMilleModule(imd);
3888 if (!mod->GeomParamsGlobal()) continue;
3889 mod->GetGeomParamsLoc(pars);
3890 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
3891 mod->SetGeomParamsGlobal(kFALSE);
3892 }
3893}
3894
3895//________________________________________________________________________________________________________
3896void AliITSAlignMille2::SetBField(Double_t b)
3897{
3898 // set Bz value
3899 if (IsZero(b,1e-5)) {
3900 fBField = 0.0;
3901 fBOn = kFALSE;
3902 fNLocal = 4;
3903 }
3904 else {
3905 fBField = b;
3906 fBOn = kTRUE;
3907 fNLocal = 5; // helices
3908 }
3909}
3910
3911//________________________________________________________________________________________________________
3912Int_t AliITSAlignMille2::ProcessUserInfo(TList* userInfo)
3913{
3914 // extract calibration information used for TrackPointArray creation from run info
3915 //
3916 if (!userInfo) { AliInfo("No UserInfo is provided"); return 0;}
3917 //
3918 TMap *cdbMap=0;
3919 TList* cdbList=0;
ef24eb3b 3920 TObjString *objStr,*objStr1,*keyStr;
3921 TString cdbStr;
6be22b3f 3922 AliCDBManager* man = AliCDBManager::Instance();
ef24eb3b 3923 man->SetCacheFlag(kFALSE);
6be22b3f 3924 //
3925 int run = userInfo->GetUniqueID();
8102b2c9 3926 if (run>0) SetRunID(run);
6be22b3f 3927 AliInfo(Form("UserInfo corresponds to run#%d",run));
3928 cdbMap = (TMap*)userInfo->FindObject("cdbMap");
ef24eb3b 3929 const TMap *curMap = man->GetStorageMap();
6be22b3f 3930 if (!cdbMap) {AliInfo("No CDB Map found in UserInfo");}
3931 else {
ef24eb3b 3932 if ((objStr=(TObjString*)cdbMap->GetValue("default"))) { // first set default CDB path
3933 if ((objStr1=(TObjString*)curMap->GetValue("default")) && objStr1->GetUniqueID()) {
3934 AliInfo(Form("OCDB default path from UserInfo: %s is overriden by user setting %s",objStr->GetName(),objStr1->GetName()));
3935 }
3936 else {
3937 cdbStr = objStr->GetString();
3938 man->UnsetDefaultStorage();
3939 if (man->GetRaw()) man->SetRaw(kFALSE);
3940 if (cdbStr.BeginsWith("raw://")) cdbStr = "raw://";
3941 AliInfo(Form("Default CDB Storage from UserInfo: %s",cdbStr.Data()));
3942 man->SetDefaultStorage( cdbStr.Data() ); // this may be overriden later by configuration file
3943 }
6be22b3f 3944 }
ef24eb3b 3945 if (man->GetRaw() && run>0) man->SetRun(run);
6be22b3f 3946 //
3947 // set specific paths relevant for alignment
3948 TIter itMap(cdbMap);
3949 while( (keyStr=(TObjString*)itMap.Next()) ) {
3950 TString keyS = keyStr->GetString();
3951 if ( keyS == "default" ) continue;
ef24eb3b 3952 //
3953 TObjString* curPath = (TObjString*)curMap->GetValue(keyStr->GetName());
3954 if (curPath && curPath->GetUniqueID()) {
3955 AliInfo(Form("Storage for %s from UserInfo\n is overriden by user setting %s",keyS.Data(),curPath->GetName()));
3956 continue;
3957 }
6be22b3f 3958 man->SetSpecificStorage( keyS.Data(), cdbMap->GetValue(keyS)->GetName() );
3959 }
3960 }
3961 //
3962 cdbList = (TList*)userInfo->FindObject("cdbList");
3963 if (!cdbList) {AliInfo("No CDB List found in UserInfo");}
3964 else {
8102b2c9 3965 // Objects used for TrackPointArray production
3966 GetPathFromUserInfo(cdbList,"GRP/Geometry/Data",fIniGeomPath ,kSameInitGeomBit);
3967 GetPathFromUserInfo(cdbList,"ITS/Align/Data" ,fIniDeltaPath,kSameInitDeltasBit);
3968 GetPathFromUserInfo(cdbList,"ITS/Calib/RespSDD",fIniSDDRespPath,kSameInitSDDRespBit);
3969 GetPathFromUserInfo(cdbList,"ITS/Calib/DriftSpeedSDD",fIniSDDVDriftPath,kSameInitSDDVDriftBit);
3970 GetPathFromUserInfo(cdbList,"ITS/Calib/MapsTimeSDD",fIniSDDCorrMapPath,kSameInitSDDCorrMapBit);
3971 GetPathFromUserInfo(cdbList,"GRP/Calib/MeanVertexSPD",fDiamondPath,kSameDiamondBit);
6be22b3f 3972 }
3973 //
1d06ac63 3974 TList *bzlst = (TList*)userInfo->FindObject("BzkGauss");
3975 if (bzlst && bzlst->At(0)) {
3976 objStr = (TObjString*)bzlst->At(0);
6be22b3f 3977 SetBField( objStr->GetString().Atof() );
8102b2c9 3978 AliInfo(Form("Magnetic field from UserInfo: %+.2e",GetBField()));
6be22b3f 3979 }
3980 return 0;
3981}
3982
8102b2c9 3983//________________________________________________________________________________________________________
b80c197e 3984Int_t AliITSAlignMille2::GetPathFromUserInfo(const TList* cdbList,const char* calib,TString& path, Int_t useBit)
8102b2c9 3985{
3986 // extract the path for specific CDB path from user info. If it is the same as already loaded, set corresponing bit
3987 TIter itList(cdbList);
3988 if (useBit>=0) ResetBit(useBit);
3989 TObjString* objStr;
3990 while( (objStr=(TObjString*)itList.Next()) )
3991 if (objStr->GetString().Contains(calib)) {
3992 TString newpath = objStr->GetString();
3993 AliInfo(Form("Found path in UserInfo: %s",newpath.Data()));
3994 if ( useBit>=0 && (fUserProvided&useBit) ) {
3995 AliInfo(Form("Will use the one provided in config: %s",path.Data()));
3996 SetBit(useBit);
3997 }
3998 else if ( useBit>=0 && (newpath == path) ) {
3999 AliInfo(Form("Path %s is the same as already loaded",path.Data()));
4000 SetBit(useBit);
4001 }
4002 else path = newpath;
4003 //
4004 return 0;
4005 }
4006 AliInfo(Form("Did not find path for %s in UserInfo",calib));
4007 path = "";
4008 return -1;
4009}
4010
6be22b3f 4011//________________________________________________________________________________________________________
4012Int_t AliITSAlignMille2::LoadSDDResponse(TString& path, AliITSresponseSDD *&resp)
4013{
8102b2c9 4014 // load SDD response
6be22b3f 4015 if (path.IsNull()) return 0;
ef24eb3b 4016 AliInfo(Form("Loading SDD response from %s",path.Data()));
6be22b3f 4017 //
abd7ef79 4018 AliCDBEntry *entry = 0;
ef24eb3b 4019 delete resp;
6be22b3f 4020 resp = 0;
8102b2c9 4021 //
6be22b3f 4022 while(1) {
4023 if (path.BeginsWith("path: ")) { // must load from OCDB
ef24eb3b 4024 entry = GetCDBEntry(path.Data());
6be22b3f 4025 if (!entry) break;
4026 resp = (AliITSresponseSDD*) entry->GetObject();
4027 entry->SetObject(NULL);
4028 entry->SetOwner(kTRUE);
8102b2c9 4029 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
ef24eb3b 4030 // delete cdbId;
4031 // delete entry;
6be22b3f 4032 break;
4033 }
4034 //
4035 if (gSystem->AccessPathName(path.Data())) break;
4036 TFile* precf = TFile::Open(path.Data());
abd7ef79 4037 if (precf->FindKey("AliITSresponseSDD")) resp = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
4038 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4039 resp = (AliITSresponseSDD*) entry->GetObject();
4040 if (resp && resp->InheritsFrom(AliITSresponseSDD::Class())) entry->SetObject(NULL);
4041 else resp = 0;
ef24eb3b 4042 entry->SetObject(NULL);
abd7ef79 4043 entry->SetOwner(kTRUE);
4044 delete entry;
4045 }
4046 //
6be22b3f 4047 precf->Close();
4048 delete precf;
4049 break;
4050 }
4051 //
4052 if (!resp) {AliError(Form("Failed to load SDD response from %s",path.Data())); return -1;}
4053 return 0;
4054}
4055
ef24eb3b 4056//________________________________________________________________________________________________________
4057Int_t AliITSAlignMille2::LoadSDDVDrift(TString& path, TObjArray *&arr)
4058{
8102b2c9 4059 // load VDrift object
ef24eb3b 4060 if (path.IsNull()) return 0;
4061 AliInfo(Form("Loading SDD VDrift from %s",path.Data()));
4062 //
4063 AliCDBEntry *entry = 0;
4064 delete arr;
4065 arr = 0;
4066 while(1) {
4067 if (path.BeginsWith("path: ")) { // must load from OCDB
4068 entry = GetCDBEntry(path.Data());
4069 if (!entry) break;
4070 arr = (TObjArray*) entry->GetObject();
4071 entry->SetObject(NULL);
4072 entry->SetOwner(kTRUE);
8102b2c9 4073 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
ef24eb3b 4074 // delete cdbId;
4075 // delete entry;
4076 break;
4077 }
4078 //
4079 if (gSystem->AccessPathName(path.Data())) break;
4080 TFile* precf = TFile::Open(path.Data());
4081 if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4082 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4083 arr = (TObjArray*) entry->GetObject();
4084 if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4085 else arr = 0;
4086 entry->SetObject(NULL);
4087 entry->SetOwner(kTRUE);
4088 delete entry;
4089 }
4090 //
4091 precf->Close();
4092 delete precf;
4093 break;
4094 }
4095 //
4096 if (!arr) {AliError(Form("Failed to load SDD vdrift from %s",path.Data())); return -1;}
4097 arr->SetOwner(kTRUE);
4098 return 0;
4099}
4100
8102b2c9 4101//________________________________________________________________________________________________________
4102Int_t AliITSAlignMille2::LoadSDDCorrMap(TString& path, AliITSCorrectSDDPoints *&map)
4103{
4104 // Load SDD correction map
4105 //
4106 if (path.IsNull()) return 0;
4107 AliInfo(Form("Loading SDD Correction Maps from %s",path.Data()));
4108 //
4109 AliCDBEntry *entry = 0;
4110 delete map;
4111 map = 0;
4112 TObjArray* arr = 0;
4113 while(1) {
4114 if (path.BeginsWith("path: ")) { // must load from OCDB
4115 entry = GetCDBEntry(path.Data());
4116 if (!entry) break;
4117 arr = (TObjArray*) entry->GetObject();
4118 entry->SetObject(NULL);
4119 entry->SetOwner(kTRUE);
4120 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
4121 // delete cdbId;
4122 // delete entry;
4123 break;
4124 }
4125 //
4126 if (gSystem->AccessPathName(path.Data())) break;
4127 TFile* precf = TFile::Open(path.Data());
4128 if (precf->FindKey("TObjArray")) arr = (TObjArray*)precf->Get("TObjArray");
4129 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4130 arr = (TObjArray*) entry->GetObject();
4131 if (arr && arr->InheritsFrom(TObjArray::Class())) entry->SetObject(NULL);
4132 else arr = 0;
4133 entry->SetObject(NULL);
4134 entry->SetOwner(kTRUE);
4135 delete entry;
4136 }
4137 //
4138 precf->Close();
4139 delete precf;
4140 break;
4141 }
4142 //
4143 if (!arr) {AliError(Form("Failed to load SDD Correction Map from %s",path.Data())); return -1;}
4144 arr->SetOwner(kTRUE);
4145 map = new AliITSCorrectSDDPoints(arr);
4146
4147 return 0;
4148}
4149
ef24eb3b 4150//________________________________________________________________________________________________________
4151Int_t AliITSAlignMille2::LoadDiamond(TString& path)
4152{
8102b2c9 4153 // load vertex constraint
ef24eb3b 4154 if (path.IsNull()) return 0;
4155 AliInfo(Form("Loading Diamond Constraint from %s",path.Data()));
4156 //
4157 AliCDBEntry *entry = 0;
4158 AliESDVertex *vtx = 0;
4159 while(1) {
4160 if (path.BeginsWith("path: ")) { // must load from OCDB
4161 entry = GetCDBEntry(path.Data());
4162 if (!entry) break;
4163 vtx = (AliESDVertex*) entry->GetObject();
4164 entry->SetObject(NULL);
4165 entry->SetOwner(kTRUE);
8102b2c9 4166 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
ef24eb3b 4167 // delete cdbId;
4168 // delete entry;
4169 break;
4170 }
4171 //
4172 if (gSystem->AccessPathName(path.Data())) break;
4173 TFile* precf = TFile::Open(path.Data());
4174 if (precf->FindKey("AliESDVertex")) vtx = (AliESDVertex*)precf->Get("AliESDVertex");
4175 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4176 vtx = (AliESDVertex*) entry->GetObject();
4177 if (vtx && vtx->InheritsFrom(AliESDVertex::Class())) entry->SetObject(NULL);
4178 else vtx = 0;
4179 entry->SetObject(NULL);
4180 entry->SetOwner(kTRUE);
4181 delete entry;
4182 }
4183 //
4184 precf->Close();
4185 delete precf;
4186 break;
4187 }
4188 //
4189 if (!vtx) {AliError(Form("Failed to load Diamond constraint from %s",path.Data())); return -1;}
4190 //
8102b2c9 4191 double vtxXYZ[3];
4192 vtx->GetXYZ(vtxXYZ);
4193 for (int i=3;i--;) vtxXYZ[i] -= fCorrDiamond[i];
4194 vtx->SetXYZ(vtxXYZ);
4195 SetVertexConstraint(vtx);
ef24eb3b 4196 AliInfo("Will use following Diamond Constraint (errors inverted):");
4197 fDiamondI.Print("");
4198 delete vtx;
4199 return 0;
4200}
4201
6be22b3f 4202//________________________________________________________________________________________________________
4203Int_t AliITSAlignMille2::LoadDeltas(TString& path, TClonesArray *&arr)
4204{
8102b2c9 4205 // load ITS geom deltas
6be22b3f 4206 if (path.IsNull()) return 0;
ef24eb3b 4207 AliInfo(Form("Loading Alignment Deltas from %s",path.Data()));
6be22b3f 4208 //
abd7ef79 4209 AliCDBEntry *entry = 0;
ef24eb3b 4210 delete arr;
6be22b3f 4211 arr = 0;
4212 while(1) {
4213 if (path.BeginsWith("path: ")) { // must load from OCDB
ef24eb3b 4214 entry = GetCDBEntry(path.Data());
6be22b3f 4215 if (!entry) break;
4216 arr = (TClonesArray*) entry->GetObject();
4217 entry->SetObject(NULL);
4218 entry->SetOwner(kTRUE);
8102b2c9 4219 // AliCDBManager::Instance()->UnloadFromCache(path); // don't want cached object, read new copy
ef24eb3b 4220 // delete cdbId;
4221 // delete entry;
6be22b3f 4222 break;
4223 }
4224 //
4225 if (gSystem->AccessPathName(path.Data())) break;
4226 TFile* precf = TFile::Open(path.Data());
abd7ef79 4227 if (precf->FindKey("ITSAlignObjs")) arr = (TClonesArray*)precf->Get("ITSAlignObjs");
4228 else if (precf->FindKey("AliCDBEntry") && (entry=(AliCDBEntry*)precf->Get("AliCDBEntry"))) {
4229 arr = (TClonesArray*) entry->GetObject();
4230 if (arr && arr->InheritsFrom(TClonesArray::Class())) entry->SetObject(NULL);
4231 else arr = 0;
ef24eb3b 4232 entry->SetObject(NULL);
abd7ef79 4233 entry->SetOwner(kTRUE);
4234 delete entry;
4235 }
6be22b3f 4236 precf->Close();
4237 delete precf;
4238 break;
4239 }
4240 //
4241 if (!arr) {AliError(Form("Failed to load Deltas from %s",path.Data())); return -1;}
ef24eb3b 4242 //
6be22b3f 4243 return 0;
4244}
4245
4246//________________________________________________________________________________________________________
abd7ef79 4247Int_t AliITSAlignMille2::CacheMatricesCurr()
6be22b3f 4248{
4249 // build arrays for the fast access to sensor matrices from their sensor ID
4250 //
4251 TGeoHMatrix mdel;
abd7ef79 4252 AliInfo("Building sensors current matrices cache");
6be22b3f 4253 //
6be22b3f 4254 fCacheMatrixCurr.Delete();
6be22b3f 4255 for (int idx=0;idx<=kMaxITSSensID;idx++) {
4256 int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
6be22b3f 4257 TGeoHMatrix *mcurr = new TGeoHMatrix();
abd7ef79 4258 AliITSAlignMille2Module::SensVolMatrix(volID, mcurr);
6be22b3f 4259 fCacheMatrixCurr.AddAtAndExpand(mcurr,idx);
4260 //
4261 }
4262 //
ef24eb3b 4263 TGeoHMatrix *mcurr = new TGeoHMatrix();
4264 fCacheMatrixCurr.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4265 //
6be22b3f 4266 fCacheMatrixCurr.SetOwner(kTRUE);
abd7ef79 4267 return 0;
4268}
4269
4270//________________________________________________________________________________________________________
4271Int_t AliITSAlignMille2::CacheMatricesOrig()
4272{
4273 // build arrays for the fast access to sensor original matrices (used for production)
4274 //
4275 TGeoHMatrix mdel;
4276 AliInfo("Building sensors original matrices cache");
4277 //
8102b2c9 4278 /*if (fIniGeomPath!=fGeometryPath)*/ if (LoadGeometry(fIniGeomPath)) {AliInfo("Failed to re-load ideal geometry");exit(1);}
4279 //
abd7ef79 4280 fCacheMatrixOrig.Delete();
ef24eb3b 4281 if (!fIniDeltaPath.IsNull()) {
4282 TClonesArray* prealSav = fPrealignment;
4283 fPrealignment = 0;
4284 if (LoadDeltas(fIniDeltaPath,fPrealignment) || ApplyToGeometry())
abd7ef79 4285 { AliError("Failed to load/apply initial deltas used to produce points"); return -1;}
ef24eb3b 4286 delete fPrealignment;
4287 fPrealignment = prealSav;
abd7ef79 4288 }
4289 //
4290 for (int idx=0;idx<=kMaxITSSensID;idx++) {
4291 int volID = AliITSAlignMille2Module::GetVolumeIDFromIndex(idx);
4292 TGeoHMatrix *morig = new TGeoHMatrix();
ef24eb3b 4293 AliITSAlignMille2Module::SensVolMatrix(volID,morig);
abd7ef79 4294 fCacheMatrixOrig.AddAtAndExpand(morig,idx);
4295 }
4296 //
8102b2c9 4297 if (fConvertPreDeltas) {
4298 // in order to convert deltas from old to new geometry we need the final matrices for all alignable objects
4299 int nmat = fGeoManager->GetNAlignable();
4300 fConvAlgMatOld.Delete();
4301 int nmatSel = 0;
4302 for (int i=0;i<nmat;i++) {
4303 TString nm = fGeoManager->GetAlignableEntry(i)->GetName();
4304 if (!nm.BeginsWith("ITS")) continue;
4305 TGeoHMatrix *mo = new TGeoHMatrix();
4306 (*mo) = *(AliGeomManager::GetMatrix(nm));
4307 fConvAlgMatOld.AddAtAndExpand(mo,nmatSel++);
4308 mo->SetTitle(nm);
4309 mo->SetName(nm);
4310 }
4311 ConvSortHierarchically(fConvAlgMatOld);
4312 }
ef24eb3b 4313 //
4314 TGeoHMatrix *mcurr = new TGeoHMatrix();
4315 fCacheMatrixOrig.AddAtAndExpand(mcurr,kVtxSensID); // special unit matrix for diamond constraint
4316 //
abd7ef79 4317 fCacheMatrixOrig.SetOwner(kTRUE);
ef24eb3b 4318
4319 fUsePreAlignment = 0;
8102b2c9 4320 LoadGeometry(fGeometryPath); // reload target geometry
ef24eb3b 4321 //
6be22b3f 4322 return 0;
4323}
4324
1d06ac63 4325//________________________________________________________________________________________________________
4326void AliITSAlignMille2::RemoveHelixFitConstraint()
4327{
4328 // suppress constraint
4329 fConstrCharge = 0;
4330 fConstrPT = fConstrPTErr = -1;
4331}
4332
6be22b3f 4333//________________________________________________________________________________________________________
4334void AliITSAlignMille2::ConstrainHelixFitPT(Int_t q,Double_t pt,Double_t pterr)
4335{
4336 // constrain q and pT of the helical fit of the track (should be set before process.track)
4337 //
4338 fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4339 fConstrPT = pt;
4340 fConstrPTErr = pterr;
8693b6b7 4341 fCurvFitWasConstrained = kTRUE;
6be22b3f 4342}
4343
4344//________________________________________________________________________________________________________
4345void AliITSAlignMille2::ConstrainHelixFitCurv(Int_t q,Double_t crv,Double_t crverr)
4346{
4347 // constrain charge and curvature of the helical fit of the track (should be set before process.track)
4348 //
4349 const double kCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
4350
4351 fConstrCharge = q==0 ? q:TMath::Sign(1,q);
4352 if (crv<0 || IsZero(crv)) {
4353 fConstrPT = -1;
4354 fConstrPTErr = -1;
8693b6b7 4355 fCurvFitWasConstrained = kFALSE;
6be22b3f 4356 }
4357 else {
8102b2c9 4358 fConstrPT = TMath::Abs(1./crv*fBField*kCQConv);
4359 fConstrPTErr = crverr>1e-10 ? TMath::Abs(fConstrPT/crv*crverr) : 0.;
8693b6b7 4360 fCurvFitWasConstrained = kTRUE;
6be22b3f 4361 }
4362}
4363
4364//________________________________________________________________________________________________________
4365TClonesArray* AliITSAlignMille2::CreateDeltas()
4366{
4367 // Create \Deltas for every explicitly or implicitly (via non-alignable volumes) varied
4368 // or prealigned module.
4369 // If the module has inded J in the hierarchy of alignable volumes (0 - the top, most
4370 // coarse level), then its Delta is expressed via MP2 \deltas (in global frame) and
4371 // prealignment \DeltaP's as:
4372 // \Delta_J = Y X Y^-1
4373 // where X = \delta_J * \DeltaP_J
4374 // Y = Prod_{K=0,J-1} \delta_K
4375 // Note that \delta_L accounts not only for its own correction but also of all non-alignable
4376 // modules in the hierarchy chain from L up to the closest alignable:
4377 // while (parent && !parent->IsAlignable()) {
4378 // \delta_L->MultiplyLeft( \delta_parent );
4379 // parent = parent->GetParent();
4380 // }
4381 //
4382 Bool_t convLoc = kFALSE;
4383 if (!GetUseGlobalDelta()) {
4384 ConvertParamsToGlobal();
4385 convLoc = kTRUE;
4386 }
4387 //
4388 AliAlignObjParams tempAlignObj;
4389 TGeoHMatrix tempMatX,tempMatY,tempMat1;
4390 //
4391 TClonesArray *array = new TClonesArray("AliAlignObjParams",10);
4392 TClonesArray &alobj = *array;
4393 int idx = 0;
4394 //
4395 TGeoManager* geoManager = AliGeomManager::GetGeometry();
4396 int nalgtot = geoManager->GetNAlignable();
4397 //
4398 for (int ialg=0;ialg<nalgtot;ialg++) { // loop over all alignable entries
4399 //
4400 const char* algname = geoManager->GetAlignableEntry(ialg)->GetName();
4401 //
4402 AliITSAlignMille2Module* md = GetMilleModuleBySymName(algname); // explicitly varied?
4403 AliITSAlignMille2Module* parent = md ? md->GetParent(): GetMilleModuleIfContained(algname);
8102b2c9 4404 if (md && parent) {
4405 TString mdName = md->GetName();
4406 TString prName = parent->GetName();
4407 // SPD Sector -> Layer parentship is fake, need special treatment
4408 if ( mdName.CountChar('/')==2 && mdName.BeginsWith("ITS/SPD") && // SPD sector
4409 prName.CountChar('/')==1 && mdName.BeginsWith("ITS/SPD") ) // SPD Layer
c82cdb65 4410 parent = parent->GetParent();//: GetMilleModuleIfContained(prName.Data());
8102b2c9 4411 }
4412 //
6be22b3f 4413 AliAlignObjParams* preob = GetPrealignedObject(algname); // was it prealigned ?
4414 //
4415 if (!preob && !md && (!parent || parent->IsAlignable())) continue; // noting to do
4416 //
4417 // create matrix X (see comment) ------------------------------------------------->>>
4418 // start from unity matrix
4419 tempMatX.Clear();
4420 if (preob) { // account prealigngment
4421 preob->GetMatrix(tempMat1);
4422 tempMatX.MultiplyLeft(&tempMat1);
4423 }
4424 //
4425 if (md) {
4426 tempAlignObj.SetTranslation( md->GetParVal(0),md->GetParVal(1),md->GetParVal(2));
4427 tempAlignObj.SetRotation( md->GetParVal(3),md->GetParVal(4),md->GetParVal(5));
4428 tempAlignObj.GetMatrix(tempMat1);
4429 tempMatX.MultiplyLeft(&tempMat1); // acount correction to varied module
4430 }
4431 //
4432 // the corrections to all non-alignable modules from current on
4433 // till first alignable should add up to its matrix
4434 while (parent && !parent->IsAlignable()) {
4435 tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4436 tempAlignObj.SetRotation( parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4437 tempAlignObj.GetMatrix(tempMat1);
4438 tempMatX.MultiplyLeft(&tempMat1); // add matrix of non-alignable module
4439 parent = parent->GetParent();
4440 }
4441 // create matrix X (see comment) ------------------------------------------------<<<
4442 //
4443 // create matrix Y (see comment) ------------------------------------------------>>>
4444 // start from unity matrix
4445 tempMatY.Clear();
4446 while ( parent ) {
4447 tempAlignObj.SetTranslation( parent->GetParVal(0),parent->GetParVal(1),parent->GetParVal(2));
4448 tempAlignObj.SetRotation( parent->GetParVal(3),parent->GetParVal(4),parent->GetParVal(5));
4449 tempAlignObj.GetMatrix(tempMat1);
4450 tempMatY.MultiplyLeft(&tempMat1);
4451 parent = parent->GetParent();
4452 }
4453 // create matrix Y (see comment) ------------------------------------------------<<<
4454 //
4455 tempMatX.MultiplyLeft(&tempMatY);
4456 tempMatX.Multiply(&tempMatY.Inverse());
4457 //
8fd71c0a 4458 if (tempMatX.IsIdentity()) continue; // do not store dummy matrices
6be22b3f 4459 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname);
4460 new(alobj[idx++]) AliAlignObjParams(algname,vid,tempMatX,kTRUE);
4461 //
4462 }
4463 //
4464 if (convLoc) ConvertParamsToLocal();
4465 //
4466 return array;
4467 //
4468}
4469
4470//_______________________________________________________________________________________
4471AliITSresponseSDD* AliITSAlignMille2::CreateSDDResponse()
4472{
4473 // create object with SDD repsonse (t0 and vdrift corrections) accounting for
4474 // eventual precalibration
4475 //
4476 // if there was a precalibration provided, copy it to new arrray
ef24eb3b 4477 AliITSresponseSDD *precal = GetSDDPrecalResp();
08517ba6 4478 if (!precal && fIniVDriftSDD) precal = GetSDDInitResp(); // InitResp is used only when IniVDrift is provided
ef24eb3b 4479 Bool_t isPreCalMult = precal&&precal->IsVDCorrMult() ? kTRUE : kFALSE;
6be22b3f 4480 AliITSresponseSDD *calibSDD = new AliITSresponseSDD();
ef24eb3b 4481 calibSDD->SetVDCorrMult(fIsSDDVDriftMult);
4482 //
4483 // copy initial values to the new object
4484 if (precal) {
4485 calibSDD->SetTimeOffset(precal->GetTimeOffset());
4486 calibSDD->SetADC2keV(precal->GetADC2keV());
4487 calibSDD->SetChargevsTime(precal->GetChargevsTime());
4488 for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) {
4489 calibSDD->SetModuleTimeZero(ind, precal->GetTimeZero(ind));
08517ba6 4490 calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kFALSE),kFALSE); // left
4491 calibSDD->SetDeltaVDrift(ind, precal->GetDeltaVDrift(ind,kTRUE ),kTRUE); // right
ef24eb3b 4492 calibSDD->SetADCtokeV(ind,precal->GetADCtokeV(ind));
4493 }
6be22b3f 4494 }
ef24eb3b 4495 else for (int ind=kSDDoffsID;ind<kSDDoffsID+kNSDDmod;ind++) calibSDD->SetModuleTimeZero(ind,0);
6be22b3f 4496 //
4497 Bool_t save = kFALSE;
4498 for (int imd=GetNModules();imd--;) {
4499 AliITSAlignMille2Module* md = GetMilleModule(imd);
4500 if (!md->IsSDD()) continue;
ef24eb3b 4501 if (md->IsFreeDOF(AliITSAlignMille2Module::kDOFT0) ||
4502 md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) ||
4503 md->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR)) save = kTRUE;
4504 //
6be22b3f 4505 for (int is=0;is<md->GetNSensitiveVolumes();is++) {
4506 int ind = md->GetSensVolIndex(is);
ef24eb3b 4507 float t0 = calibSDD->GetTimeZero(ind) + md->GetParVal(AliITSAlignMille2Module::kDOFT0);
4508 double dvL = md->GetParVal(AliITSAlignMille2Module::kDOFDVL);
4509 double dvR = md->GetParVal(AliITSAlignMille2Module::kDOFDVR);
4510 if (!calibSDD->IsVDCorrMult()) { // save as additive correction
4511 dvL *= 1e4;
4512 dvR *= 1e4;
4513 //
4514 double conv = 1;
4515 if (isPreCalMult) conv = 6.4; // convert multiplicative precal correction to additive
4516 dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4517 dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4518 }
4519 else { // save as multipicative correction
4520 double conv = 1;
4521 if (!isPreCalMult) conv = 1./6.4; // convert additive precal correction to multiplicative
4522 dvL += calibSDD->GetDeltaVDrift(ind,kFALSE)*conv;
4523 dvR += calibSDD->GetDeltaVDrift(ind,kTRUE)*conv;
4524 }
6be22b3f 4525 //
4526 calibSDD->SetModuleTimeZero(ind, t0);
ef24eb3b 4527 calibSDD->SetDeltaVDrift(ind, dvL, kFALSE); // left side correction
4528 calibSDD->SetDeltaVDrift(ind, dvR, kTRUE); // right side correction
6be22b3f 4529 }
4530 }
4531 //
4532 if (!save) {
4533 AliInfo("No free parameters for SDD calibration, nothing to save");
4534 delete calibSDD;
4535 calibSDD = 0;
4536 }
4537 //
4538 return calibSDD;
4539}
24391cd5 4540
ef24eb3b 4541//_______________________________________________________________________________________
4542Int_t AliITSAlignMille2::ReloadInitCalib(TList *userInfo)
4543{
4544 // Use provided UserInfo to
4545 // load the initial calib parameters (geometry, SDD response...)
4546 // Can be used if set of data was processed with different calibration
4547 //
4548 if (!userInfo) {
4549 AliInfo("Reloading of the Calibration parameters was called with empty userInfo");
4550 return 1;
4551 }
4552 if (ProcessUserInfo(userInfo)) {
4553 AliInfo("Error in processing user info");
4554 userInfo->Print();
4555 exit(1);
4556 }
4557 return ReloadInitCalib();
4558}
4559
4560//_______________________________________________________________________________________
4561Int_t AliITSAlignMille2::ReloadInitCalib()
4562{
4563 // Load the initial calib parameters (geometry, SDD response...)
4564 // Can be used if set of data was processed with different calibration
4565 //
4566 // 1st cache original matrices
8102b2c9 4567 if (!(TestBit(kSameInitDeltasBit) && TestBit(kSameInitGeomBit))) { // need to reload geometry
4568 //
ef24eb3b 4569 if (CacheMatricesOrig()) {
4570 AliInfo("Failed to cache new initial geometry");
4571 exit(1);
4572 }
8102b2c9 4573 // RS : commented because we don't need to reload prealignment deltas, they are already loaded
4574 // then reload the prealignment geometry
4575 // if (LoadDeltas(fPreDeltaPath,fPrealignment)) {
4576 // AliInfo(Form("Failed to reload the prealigned geometry %s",fPreDeltaPath.Data()));
4577 // exit(1);
4578 // }
ef24eb3b 4579 //
4580 if (fPrealignment && ApplyToGeometry()) {
4581 AliInfo(Form("Failed re-apply prealigned geometry %s",fPreDeltaPath.Data()));
4582 exit(1);
4583 }
4584 //
4585 // usually no need to re-cache the prealignment geometry, it was not changed
4586 if (fCacheMatrixCurr.GetEntriesFast() != fCacheMatrixOrig.GetEntriesFast()) {
4587 // CacheMatricesCurr();
4588 AliInfo(Form("Failed to cache the prealigned geometry %s",fPreDeltaPath.Data()));
4589 exit(1);
4590 }
4591 }
4592 else ResetBit(kSameInitDeltasBit);
4593 //
4594 // reload initial SDD response
4595 if (!TestBit(kSameInitSDDRespBit)) {
4596 if (LoadSDDResponse(fIniSDDRespPath, fIniRespSDD) ) {
4597 AliInfo(Form("Failed to load new SDD response %s",fIniSDDRespPath.Data()));
4598 exit(1);
4599 }
4600 }
4601 else ResetBit(kSameInitSDDRespBit);
4602 //
4603 // reload initial SDD vdrift
4604 if (!TestBit(kSameInitSDDVDriftBit)) {
4605 if (LoadSDDVDrift(fIniSDDVDriftPath, fIniVDriftSDD) ) {
4606 AliInfo(Form("Failed to load new SDD VDrift %s",fIniSDDVDriftPath.Data()));
4607 exit(1);
4608 }
4609 }
4610 else ResetBit(kSameInitSDDRespBit);
4611 //
8102b2c9 4612 // reload SDD corr.map
4613 if (!TestBit(kSameInitSDDCorrMapBit)) {
4614 if (LoadSDDCorrMap(fIniSDDCorrMapPath, fIniCorrMapSDD) ) {
4615 AliInfo(Form("Failed to load new SDD Correction Map %s",fIniSDDCorrMapPath.Data()));
4616 exit(1);
4617 }
4618 }
4619 else ResetBit(kSameInitSDDRespBit);
4620 //
ef24eb3b 4621 // reload diamond info
4622 if (!TestBit(kSameDiamondBit)) {
4623 if (LoadDiamond(fDiamondPath) ) {
4624 AliInfo(Form("Failed to load new Diamond constraint %s",fDiamondPath.Data()));
4625 exit(1);
4626 }
4627 }
4628 else ResetBit(kSameInitSDDRespBit);
4629 //
ef24eb3b 4630 return 0;
4631}
4632
4633//_______________________________________________________________________________________
4634void AliITSAlignMille2::JacobianPosGloLoc(int locid,double* jacobian)
4635{
4636 // calculate the locid row of the jacobian for transformation of the local coordinate to global at current point
4637 TGeoHMatrix* mat = GetSensorCurrMatrixSID(fCurrentSensID);
4638 const Double_t dpar = 1e-2;
4639 double sav = fMeasLoc[locid];
4640 fMeasLoc[locid] += dpar;
4641 mat->LocalToMaster(fMeasLoc,jacobian);
4642 fMeasLoc[locid] = sav; // recover original value
4643 for (int i=3;i--;) jacobian[i] = (jacobian[i]-fMeasGlo[i])/dpar; // the transformation is linear!!!
4644}
4645
4646//_______________________________________________________________________________________
4647void AliITSAlignMille2::TieSDDVDriftsLR(AliITSAlignMille2Module* mod)
4648{
4649 // impose equality of Left/Right sides VDrift correction for SDD
4650 ResetLocalEquation();
4651 if ( (mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVL) + mod->IsFreeDOF(AliITSAlignMille2Module::kDOFDVR))==1) {
4652 AliError("Left/Right VDrift equality is requested for SDD module with only one side VDrift free");
4653 mod->Print();
4654 return;
4655 }
8102b2c9 4656 if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVL), 1.);
4657 if (mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR)>=0) SetGlobalDerivative(mod->GetParOffset(AliITSAlignMille2Module::kDOFDVR), -1.);
ef24eb3b 4658 AddConstraint(fGlobalDerivatives, 0, 1e-12);
4659 //
4660}
4661
4662//_______________________________________________________________________________________
4663void AliITSAlignMille2::ProcessSDDPointInfo(const AliTrackPoint* pnt,Int_t sID, Int_t pntID)
4664{
4665 // extract the drift information from SDD track point
4666 //
4667 fDriftTime0[pntID] = fIniRespSDD ? fIniRespSDD->GetTimeZero(sID) : 0.;
4668 double tdif = pnt->GetDriftTime() - fDriftTime0[pntID];
4669 if (tdif<0) tdif = 1;
4670 //
4671 // VDrift extraction
08517ba6 4672 double vdrift=0,vdrift0=0;
ef24eb3b 4673 Bool_t sddSide = kFALSE;
4674 int sID0 = 2*(sID-kSDDoffsID);
4675 double zanode = -999;
4676 //
4677 if (fIniVDriftSDD) { // SDD VDrift object is provided, use the vdrift from it
4678 AliITSDriftSpeedArraySDD* drarr;
4679 double vdR,vdL,xlR,xlL;
4680 // sometimes xlocal on right side is negative due to the wrong calibration, need to test both hypothesis
4681 double xlabs = TMath::Abs(fMeasLoc[kX]);
4682 drarr = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0); // left side, xloc>0
4683 zanode = fSegmentationSDD->GetAnodeFromLocal(xlabs,fMeasLoc[kZ]);
4684 vdL = drarr->GetDriftSpeed(0, zanode);
4685 if (fIniRespSDD) {
4686 double corr = fIniRespSDD->GetDeltaVDrift(sID, kFALSE);
4687 if (fIniRespSDD->IsVDCorrMult()) vdL *= (1+corr);
4688 else vdL += corr;
4689 }
4690 xlL = (fSegmentationSDD->Dx() - vdL*tdif)*1e-4;
4691 //
4692 drarr = (AliITSDriftSpeedArraySDD*)fIniVDriftSDD->At(sID0+1); // right side, xloc<0
4693 zanode = fSegmentationSDD->GetAnodeFromLocal(-xlabs,fMeasLoc[kZ]) - 256;
4694 vdR = drarr->GetDriftSpeed(0, zanode);
4695 if (fIniRespSDD) {
4696 double corr = fIniRespSDD->GetDeltaVDrift(sID, kTRUE);
4697 if (fIniRespSDD->IsVDCorrMult()) vdR *= (1+corr);
4698 else vdR += corr;
4699 }
4700 xlR = -(fSegmentationSDD->Dx() - vdR*tdif)*1e-4;
4701 //
4702 if (TMath::Abs(xlL-fMeasLoc[kX])<TMath::Abs(xlR-fMeasLoc[kX])) {
4703 sddSide = 0; // left side
4704 vdrift = vdL*1e-4;
4705 }
4706 else { // right side
4707 sddSide = 1;
4708 vdrift = vdR*1e-4;
4709 }
4710 //
4711 }
4712 else { // try to determine the vdrift from the xloc
4713 vdrift = (fSegmentationSDD->Dx()*1e-4 - TMath::Abs(fMeasLoc[kX]))/tdif;
4714 sddSide = fMeasLoc[kX]<0; // 0 = left (xloc>0) ; 1 = right (xloc<1)
4715 }
4716 //
4717 if (fPreVDriftSDD) { // use imposed vdrift as a starting point
4718 zanode = fSegmentationSDD->GetAnodeFromLocal(0.5-sddSide,fMeasLoc[kZ]);
4719 if (sddSide) zanode -= 256;
4720 vdrift = ((AliITSDriftSpeedArraySDD*)fPreVDriftSDD->At(sID0+sddSide))->GetDriftSpeed(0, zanode)*1e-4;
4721 }
4722 //
4723 if (vdrift<0) vdrift = 0;
08517ba6 4724 vdrift0 = vdrift;
ef24eb3b 4725 // at this point we have vdrift and t0 used to create the original point.
4726 // see if precalibration was provided
4727 if (fPreRespSDD) {
4728 float t0Upd = fPreRespSDD->GetTimeZero(sID);
4729 double corr = fPreRespSDD->GetDeltaVDrift(sID, sddSide);
4730 if (fPreRespSDD->IsVDCorrMult()) vdrift *= 1+corr; // right side (xloc<0) may have different correction
4731 else vdrift += corr*1e-4;
08517ba6 4732 //
4733 // if IniRespSDD was used, it should be subtracted back, since it is accounted in the PreResp
921c1cc0 4734 if (fIniVDriftSDD&&fIniRespSDD && (fPreVDriftSDD==0)) {
08517ba6 4735 double corr1 = fIniRespSDD->GetDeltaVDrift(sID, sddSide);
4736 if (fIniRespSDD->IsVDCorrMult()) vdrift *= (1-corr1);
4737 else vdrift -= corr1*1e-4;
4738 }
ef24eb3b 4739 tdif = pnt->GetDriftTime() - t0Upd;
4740 // correct Xlocal
4741 fMeasLoc[0] = fSegmentationSDD->Dx()*1e-4 - vdrift*tdif;
4742 if (sddSide) fMeasLoc[0] = -fMeasLoc[0];
4743 fDriftTime0[pntID] = t0Upd;
4744 }
8102b2c9 4745 //
4746 if (fPreCorrMapSDD) { // apply correction map
4747 fMeasLoc[0] += fPreCorrMapSDD->GetCorrection(sID,fMeasLoc[2],fMeasLoc[0]);
4748 }
4749
ef24eb3b 4750 // TEMPORARY CORRECTION (if provided) --------------<<<
08517ba6 4751 fDriftSpeed[pntID] = sddSide ? -vdrift : vdrift;
4752 fDriftSpeed0[pntID] = sddSide ? -vdrift0 : vdrift0;
ef24eb3b 4753 //
4754 // printf("#%d: t:%+e x:%+e v:%+e: side:%d\n",pntID,fDriftTime0[pntID],fMeasLoc[0],fDriftSpeed[pntID],sddSide);
4755}
4756
4757//_______________________________________________________________________________________
4758AliITSAlignMille2Module* AliITSAlignMille2::CreateVertexModule()
4759{
4760 // creates dummy module for vertex constraint
4761 TGeoHMatrix mt;
4762 AliITSAlignMille2Module* mod = new AliITSAlignMille2Module(kVtxSensID,kVtxSensVID,"VTX",&mt,0,0);
4763 fMilleModule.AddAtAndExpand(mod,fNModules);
4764 mod->SetGeomParamsGlobal(fUseGlobalDelta);
4765 fDiamondModID = fNModules;
4766 mod->SetUniqueID(fNModules++);
4767 mod->SetNotInConf(kTRUE);
4768 return mod;
4769 //
4770}
4771
4772//_______________________________________________________________________________________
4773AliCDBEntry* AliITSAlignMille2::GetCDBEntry(const char* path)
4774{
4775 // return object from the OCDB
4776 AliCDBEntry *entry = 0;
4777 AliInfo(Form("Loading object %s",path));
4778 AliCDBManager* man = AliCDBManager::Instance();
4779 AliCDBId* cdbId = AliCDBId::MakeFromString(path);
4780 if (!cdbId) {
4781 AliError("Failed to create cdbId");
4782 return 0;
4783 }
4784 //
4785 AliCDBStorage* stor = man->GetDefaultStorage();
4786 if (!stor && !man->GetRaw()) man->SetDefaultStorage("raw://");
8102b2c9 4787 if (man->GetRaw()) man->SetRun(fRunID>0 ? fRunID : cdbId->GetFirstRun());
ef24eb3b 4788 if (stor) {
4789 TString tp = stor->GetType();
4790 if (tp.Contains("alien",TString::kIgnoreCase) && !gGrid) TGrid::Connect("alien:");
4791 }
8102b2c9 4792 entry = man->Get(cdbId->GetPath(),cdbId->GetFirstRun(),cdbId->GetVersion(),cdbId->GetSubVersion());
4793 // entry = man->Get( *cdbId );
ef24eb3b 4794 man->ClearCache();
4795 //
4796 delete cdbId;
4797 return entry;
4798 //
4799}
4800
8102b2c9 4801//_______________________________________________________________________________________
4802void AliITSAlignMille2::SetVertexConstraint(const AliESDVertex* vtx)
4803{
4804 // set vertex for constraint
4805 if (!vtx) return;
4806 //
4807 double cmat[6];
4808 float cmatF[6];
4809 vtx->GetCovMatrix(cmat);
4810 AliITSAlignMille2Module* diamMod = GetMilleModuleByVID(kVtxSensVID);
4811 if (diamMod) {
4812 cmat[0] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaXFactor();
4813 cmat[2] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaYFactor();
4814 cmat[5] *= diamMod->GetSigmaZFactor()*diamMod->GetSigmaZFactor();
4815 cmat[1] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaYFactor();
4816 cmat[3] *= diamMod->GetSigmaXFactor()*diamMod->GetSigmaZFactor();
4817 cmat[4] *= diamMod->GetSigmaYFactor()*diamMod->GetSigmaZFactor();
4818 }
4819 cmatF[0] = cmat[0]; // xx
4820 cmatF[1] = cmat[1]; // xy
4821 cmatF[2] = cmat[3]; // xz
4822 cmatF[3] = cmat[2]; // yy
4823 cmatF[4] = cmat[4]; // yz
4824 cmatF[5] = cmat[5]; // zz
4825
4826 fDiamond.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4827 //
4828 Double_t t0 = cmat[2]*cmat[5] - cmat[4]*cmat[4];
4829 Double_t t1 = cmat[1]*cmat[5] - cmat[3]*cmat[4];
4830 Double_t t2 = cmat[1]*cmat[4] - cmat[2]*cmat[3];
4831 Double_t det = cmat[0]*t0 - cmat[1]*t1 + cmat[3]*t2;
4832 if (TMath::Abs(det)<1e-36) {
4833 vtx->Print();
4834 AliFatal("Vertex constraint cov.matrix is singular");
4835 }
4836 cmatF[0] = t0/det;
4837 cmatF[1] = -t1/det;
4838 cmatF[2] = t2/det;
4839 cmatF[3] = (cmat[0]*cmat[5] - cmat[3]*cmat[3])/det;
4840 cmatF[4] = (cmat[1]*cmat[3] - cmat[0]*cmat[4])/det;
4841 cmatF[5] = (cmat[0]*cmat[2] - cmat[1]*cmat[1])/det;
4842 fDiamondI.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
4843 fVertexSet = kTRUE;
4844 //
4845}
4846
4847//_______________________________________________________________________________________
4848void AliITSAlignMille2::ConvertDeltas()
4849{
4850 // convert prealignment deltas from old geometry to new one
4851 // NOTE: the target geometry must be loaded at time this method is called
4852 //
4853 // NOTE: This method can be ONLY used when as a prealignment deltas those used for the production
4854 // of trackpoints (e.g. extracted from the UserInfo).
4855 // The prealignment deltas provided by user via config file must be already converted to target geometry:
4856 // this can be done externally using the macro ConvertDeltas.C
4857 //
4858 // delta_j_new = delta_j_old * Xj_old * Xj_new^-1
4859 // where X = Prod{delta_i,i=j-1:0} M_j
4860 // with j - the level of the alignable volume in the hierarchy, M - corresponding ideal matrix
4861 // Note that delta_j * Xj is equal to final (misaligned) matrix of corresponding geometry, G_j.
4862 // Since this method is used ONLY in the case where the prealignment deltas are equal to production deltas,
4863 // we have already loaded G_j_old in the fConvAlgMatOld (filled in the CacheMatricesOrig)
4864 // Hence, delta_j_new = G_j_old * Xj_new^-1
4865 //
4866 AliInfo("Converting deltas from initial to target geometry");
4867 int nMatOld = fConvAlgMatOld.GetEntriesFast(); // number of alignable matrices
4868 TClonesArray* deltArrNew = new TClonesArray("AliAlignObjParams",10);
4869 //
4870 TGeoHMatrix dmPar;
4871 int nDelNew = 0;
4872 //
4873 for (int im=0;im<nMatOld;im++) {
4874 TGeoHMatrix* mtGjold = (TGeoHMatrix*)fConvAlgMatOld[im];
4875 TString algname = mtGjold->GetTitle();
4876 UShort_t vID = AliITSAlignMille2Module::GetVolumeIDFromSymname(algname.Data());
4877 //
4878 // build X_new >>>
4879 TGeoHMatrix* parent = mtGjold;
4880 TGeoHMatrix xNew;
4881 int parID;
4882 while ( (parID=parent->GetUniqueID()-1)>=0 ) {
4883 parent = (TGeoHMatrix*)fConvAlgMatOld[parID];
4884 AliAlignObjParams* deltaPar = ConvFindDelta(deltArrNew,parent->GetTitle());
4885 if (deltaPar) deltaPar->GetMatrix(dmPar); xNew *= dmPar;
4886 }
4887 AliGeomManager::GetOrigGlobalMatrix(algname,dmPar); // ideal matrix of new geometry
4888 xNew *= dmPar;
4889 // build X_new <<<
4890 //
4891 dmPar = *mtGjold;
4892 dmPar *= xNew.Inverse();
4893 new((*deltArrNew)[nDelNew++]) AliAlignObjParams(algname.Data(),vID,dmPar,kTRUE);
4894 //
4895 }
4896 delete fPrealignment;
4897 fPrealignment = deltArrNew;
4898 //
4899 // we don't need anymore old matrices
4900 fConvAlgMatOld.Delete();
4901 //
4902}
4903
4904//_______________________________________________________________________________________
4905void AliITSAlignMille2::ConvSortHierarchically(TObjArray& matArr)
4906{
4907 // Used only for the deltas conversion from one geometry to another
4908 // Sort the matrices according to hiearachy (coarse -> fine)
4909 //
4910 int nmat = matArr.GetEntriesFast();
4911 //
4912 for (int i=0;i<nmat;i++) {
4913 for (int j=i+1;j<nmat;j++) {
4914 TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4915 TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4916 if (ConvIsJParentOfI(matI,matJ)) { // swap
4917 matArr[i] = matJ;
4918 matArr[j] = matI;
4919 }
4920 }
4921 }
4922 //
4923 // set direct parent id's in the UniqueID's
4924 for (int i=nmat;i--;) {
4925 TGeoHMatrix* matI = (TGeoHMatrix*) matArr[i];
4926 matI->SetUniqueID(0);
4927 for (int j=i;j--;) {
4928 TGeoHMatrix* matJ = (TGeoHMatrix*) matArr[j];
4929 if (ConvIsJParentOfI(matI,matJ)) { matI->SetUniqueID(j+1); break; }
4930 }
4931 }
4932}
4933
4934//_______________________________________________________________________________________
4935Bool_t AliITSAlignMille2::ConvIsJParentOfI(const TGeoHMatrix* matI,const TGeoHMatrix* matJ) const
4936{
4937 // Used only for the deltas conversion from one geometry to another
4938 // True if matJ is higher in hierarchy than
4939 //
4940 TString nmI = matI->GetTitle();
4941 TString nmJ = matJ->GetTitle();
4942 //
4943 int nlrI = nmI.CountChar('/');
4944 int nlrJ = nmJ.CountChar('/');
4945 if (nlrJ>=nlrI) return kFALSE;
4946 //
4947 // special case of SPD sectors
4948 if (nmI.BeginsWith("ITS/SPD1") && nmJ.BeginsWith("ITS/SPD0") && nlrJ==2) nmJ.ReplaceAll("SPD0","SPD1");
4949 return (nmI.BeginsWith(nmJ)) ? kTRUE:kFALSE;
4950 //
4951}
4952
4953//_______________________________________________________________________________________
4954AliAlignObjParams* AliITSAlignMille2::ConvFindDelta(const TClonesArray* arrDelta,const TString& algname) const
4955{
4956 // find the delta for given module
4957 if (!arrDelta) return 0;
4958 AliAlignObjParams* delta = 0;
4959 int nDeltas = arrDelta->GetEntries();
4960 for (int id=0;id<nDeltas;id++) {
4961 delta = (AliAlignObjParams*)arrDelta->At(id);
4962 if (algname==delta->GetSymName()) break;
4963 delta = 0;
4964 }
4965 return delta;
4966}