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