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