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