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