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