]> 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
Corrected name of the SPD fast-or minimum bias trigger input (Annalisa)
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille2.cxx
CommitLineData
7b85e477 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
057304f3 16/* $Id$ */
7b85e477 17
18//-----------------------------------------------------------------------------
6526a72c 19//
20// Interface to AliMillePede2 alignment class for the ALICE ITS detector
21//
22// ITS specific alignment class which interface to AliMillepede.
23// For each track ProcessTrack calculates the local and global derivatives
24// at each hit and fill the corresponding local equations. Provide methods for
25// fixing or constraining detection elements for best results.
26//
27// author M. Lunardon (thanks to J. Castillo), ruben.shahoyan@cern.ch
7b85e477 28//-----------------------------------------------------------------------------
29
30#include <TF1.h>
31#include <TFile.h>
32#include <TClonesArray.h>
33#include <TGraph.h>
34#include <TGeoMatrix.h>
35#include <TMath.h>
36#include <TGraphErrors.h>
37#include <TVirtualFitter.h>
6526a72c 38#include <TGeoManager.h>
7b85e477 39
40#include "AliITSAlignMille2.h"
41#include "AliITSgeomTGeo.h"
42#include "AliGeomManager.h"
43#include "AliMillePede2.h"
44#include "AliTrackPointArray.h"
45#include "AliAlignObjParams.h"
46#include "AliLog.h"
47#include "TSystem.h" // come si fa?
48#include "AliTrackFitterRieman.h"
49
6526a72c 50
7b85e477 51ClassImp(AliITSAlignMille2)
6526a72c 52
53
54//========================================================================================================
7b85e477 55
56AliITSAlignMille2* AliITSAlignMille2::fgInstance = 0;
57Int_t AliITSAlignMille2::fgInstanceID = 0;
58
6526a72c 59//________________________________________________________________________________________________________
60AliITSAlignMille2::AliITSAlignMille2(const Char_t *configFilename )
7b85e477 61: TObject(),
62 fMillepede(0),
63 fStartFac(16.),
64 fResCutInitial(100.),
65 fResCut(100.),
66 fNGlobal(0),
67 fNLocal(4),
68 fNStdDev(3),
69 fIsMilleInit(kFALSE),
6526a72c 70 fAllowPseudoParents(kFALSE),
71 //
7b85e477 72 fCurrentModule(0),
73 fTrack(0),
6526a72c 74 fTrackBuff(0),
7b85e477 75 fCluster(),
76 fGlobalDerivatives(0),
6526a72c 77 //
7b85e477 78 fMinNPtsPerTrack(3),
79 fInitTrackParamsMeth(1),
80 fTotBadLocEqPoints(0),
81 fRieman(0),
82 //
6526a72c 83 fConstraints(0),
84 //
7b85e477 85 fUseGlobalDelta(kFALSE),
86 fRequirePoints(kFALSE),
87 fTempExcludedModule(-1),
88 //
89 fGeometryFileName("geometry.root"),
90 fPreAlignmentFileName(""),
6526a72c 91 fConstrRefFileName(""),
7b85e477 92 fGeoManager(0),
93 fIsConfigured(kFALSE),
94 fPreAlignQF(0),
95//
6526a72c 96 fCorrectSDD(0),
97 fInitialRecSDD(0),
7b85e477 98 fPrealignment(0),
6526a72c 99 fConstrRef(0),
7b85e477 100 fMilleModule(2),
101 fSuperModule(2),
102 fNModules(0),
103 fNSuperModules(0),
104 fUsePreAlignment(kFALSE),
7b85e477 105 fBOn(kFALSE),
106 fBField(0.0),
6526a72c 107 fBug(0),
108 fMilleVersion(2)
7b85e477 109{
110 /// main constructor that takes input from configuration file
7b85e477 111 for (int i=3;i--;) fSigmaFactor[i] = 1.0;
112 //
113 // new RS
114 for (Int_t i=0; i<6; i++) {
115 fNReqLayUp[i]=0;
116 fNReqLayDown[i]=0;
117 fNReqLay[i]=0;
118 }
119 for (Int_t i=0; i<3; i++) {
120 fNReqDetUp[i]=0;
121 fNReqDetDown[i]=0;
122 fNReqDet[i]=0;
123 }
124 //
125 Int_t lc=LoadConfig(configFilename);
126 if (lc) {
6526a72c 127 AliError(Form("Error %d loading configuration from %s",lc,configFilename));
128 exit(1);
7b85e477 129 }
130 //
6526a72c 131 fMillepede = new AliMillePede2();
7b85e477 132 fgInstance = this;
133 fgInstanceID++;
134 //
135}
136
6526a72c 137//________________________________________________________________________________________________________
138AliITSAlignMille2::~AliITSAlignMille2()
139{
7b85e477 140 /// Destructor
6526a72c 141 if (fMillepede) delete fMillepede; fMillepede = 0;
142 if (fGlobalDerivatives) delete[] fGlobalDerivatives; fGlobalDerivatives = 0;
143 if (fRieman) delete fRieman; fRieman = 0;
144 if (fPrealignment) delete fPrealignment; fPrealignment = 0;
145 if (fConstrRef) delete fConstrRef; fConstrRef = 0;
146 if (fCorrectSDD) delete fCorrectSDD; fCorrectSDD = 0;
147 if (fInitialRecSDD) delete fInitialRecSDD; fInitialRecSDD = 0;
148 fTrackBuff.Delete();
149 fConstraints.Delete();
7b85e477 150 fMilleModule.Delete();
151 fSuperModule.Delete();
152 if (--fgInstanceID==0) fgInstance = 0;
153}
154
155///////////////////////////////////////////////////////////////////////
6526a72c 156TObjArray* AliITSAlignMille2::GetConfigRecord(FILE* stream, TString& recTitle, TString& recOpt, Bool_t rew)
157{
158 TString record;
159 static TObjArray* recElems = 0;
160 if (recElems) {delete recElems; recElems = 0;}
161 //
162 TString keyws = recTitle;
163 if (!keyws.IsNull()) {
164 keyws.ToUpper();
165 // keyws += " ";
166 }
167 while (record.Gets(stream)) {
168 int cmt=record.Index("#");
169 if (cmt>=0) record.Remove(cmt); // skip comment
170 record.ReplaceAll("\t"," ");
171 record.ReplaceAll("\r"," ");
172 record.Remove(TString::kBoth,' ');
173 if (record.IsNull()) continue; // nothing to decode
174 if (!keyws.IsNull() && !record.BeginsWith(keyws.Data())) continue; // specific record was requested
175 //
176 recElems = record.Tokenize(" ");
177 recTitle = recElems->At(0)->GetName();
178 recTitle.ToUpper();
179 recOpt = recElems->GetLast()>0 ? recElems->At(1)->GetName() : "";
180 break;
181 }
182 if (rew || !recElems) rewind(stream);
183 return recElems;
184}
7b85e477 185
6526a72c 186//________________________________________________________________________________________________________
187Int_t AliITSAlignMille2::LoadConfig(const Char_t *cfile)
188{ /// return 0 if success
7b85e477 189 /// 1 if error in module index or voluid
6526a72c 190 //
7b85e477 191 FILE *pfc=fopen(cfile,"r");
192 if (!pfc) return -1;
6526a72c 193 //
194 TString record,recTitle,recOpt,recExt;
195 Int_t nrecElems,irec;
196 TObjArray *recArr=0;
197 //
198 fNModules = 0;
199 Bool_t stopped = kFALSE;
200 //
201 while(1) {
7b85e477 202 //
6526a72c 203 // ============= 1: we read some obligatory records in predefined order ================
204 //
205 recTitle = "GEOMETRY_FILE";
206 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
207 (fGeometryFileName=recOpt).IsNull() ||
208 gSystem->AccessPathName(recOpt.Data()) ||
209 InitGeometry() )
210 { AliError("Failed to find/load Geometry"); stopped = kTRUE; break;}
7b85e477 211 //
6526a72c 212 recTitle = "SUPERMODULE_FILE";
213 if ( !GetConfigRecord(pfc,recTitle,recOpt,1) ||
214 recOpt.IsNull() ||
215 gSystem->AccessPathName(recOpt.Data()) ||
216 LoadSuperModuleFile(recOpt.Data()))
217 { AliError("Failed to find/load SuperModules"); stopped = kTRUE; break;}
7b85e477 218 //
6526a72c 219 recTitle = "CONSTRAINTS_REFERENCE_FILE"; // LOCAL_CONSTRAINTS are defined wrt these deltas
220 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
221 if (recOpt.IsNull() || recOpt=="IDEAL") SetConstraintWrtRef( "IDEAL" );
222 else if (gSystem->AccessPathName(recOpt.Data()) || SetConstraintWrtRef(recOpt.Data()) )
223 { AliError("Failed to load reference deltas for local constraints"); stopped = kTRUE; break;}
7b85e477 224 }
6526a72c 225 //
226 recTitle = "PREALIGNMENT_FILE";
227 if ( GetConfigRecord(pfc,recTitle,recOpt,1) )
228 if ( (fPreAlignmentFileName=recOpt).IsNull() ||
229 gSystem->AccessPathName(recOpt.Data()) ||
230 ApplyToGeometry())
231 { AliError(Form("Failed to load Prealignment file %s",recOpt.Data())); stopped = kTRUE; break;}
7b85e477 232 //
6526a72c 233 recTitle = "PRECALIBSDD_FILE";
234 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
235 if ( recOpt.IsNull() || gSystem->AccessPathName(recOpt.Data()) ) {stopped = kTRUE; break;}
236 AliInfo(Form("Using %s for SDD precalibration",recOpt.Data()));
237 TFile* precfi = TFile::Open(recOpt.Data());
238 if (!precfi->IsOpen()) {stopped = kTRUE; break;}
239 fCorrectSDD = (AliITSresponseSDD*)precfi->Get("AliITSresponseSDD");
240 precfi->Close();
241 delete precfi;
242 if (!fCorrectSDD) {AliError("Precalibration SDD object is not found"); stopped = kTRUE; break;}
7b85e477 243 }
244 //
6526a72c 245 recTitle = "INITCALBSDD_FILE";
246 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) {
247 if ( recOpt.IsNull() || gSystem->AccessPathName(recOpt.Data()) ) {stopped = kTRUE; break;}
248 AliInfo(Form("Using %s as SDD calibration used in TrackPoints",recOpt.Data()));
249 TFile* precf = TFile::Open(recOpt.Data());
250 if (!precf->IsOpen()) {stopped = kTRUE; break;}
251 fInitialRecSDD = (AliITSresponseSDD*)precf->Get("AliITSresponseSDD");
252 precf->Close();
253 delete precf;
254 if (!fInitialRecSDD) {AliError("Initial Calibration SDD object is not found"); stopped = kTRUE; break;}
7b85e477 255 }
256 //
6526a72c 257 recTitle = "SET_GLOBAL_DELTAS";
258 if ( GetConfigRecord(pfc,recTitle,recOpt,1) ) SetUseGlobalDelta(kTRUE);
7b85e477 259 //
6526a72c 260 // =========== 2: see if there are local gaussian constraints defined =====================
261 // Note that they should be loaded before the modules declaration
7b85e477 262 //
6526a72c 263 recTitle = "CONSTRAINT_LOCAL";
264 while( (recArr=GetConfigRecord(pfc,recTitle,recOpt,0)) ) {
265 nrecElems = recArr->GetLast()+1;
266 if (recOpt.IsFloat()) {stopped = kTRUE; break;} // wrong name
267 if (GetConstraint(recOpt.Data())) {
268 AliError(Form("Existing constraint %s repeated",recOpt.Data()));
269 stopped = kTRUE; break;
270 }
271 recExt = recArr->At(2)->GetName();
272 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
273 double val = recExt.Atof();
274 recExt = recArr->At(3)->GetName();
275 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
276 double err = recExt.Atof();
277 int nwgh = nrecElems - 4;
278 double *wgh = new double[nwgh];
279 for (nwgh=0,irec=4;irec<nrecElems;irec++) {
280 recExt = recArr->At(irec)->GetName();
281 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
282 wgh[nwgh++] = recExt.Atof();
283 }
284 if (stopped) {delete[] wgh; break;}
285 //
286 ConstrainLocal(recOpt.Data(),wgh,nwgh,val,err);
287 delete[] wgh;
288 //
289 } // end while for loop over local constraints
290 if (stopped) break;
7b85e477 291 //
6526a72c 292 // =========== 3: now read modules to align ===================================
7b85e477 293 //
6526a72c 294 rewind(pfc);
295 while( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0)) ) {
296 if (!(recTitle=="MODULE_VOLUID" || recTitle=="MODULE_INDEX")) continue;
297 // Expected format: MODULE id tolX tolY tolZ tolPsi tolTh tolPhi [[sigX sigY sigZ] extra params]
298 // where tol* is the tolerance (sigma) for given DOF. 0 means fixed
299 // sig* is the scaling parameters for the errors of the clusters of this module
300 // extra params are defined for specific modules, e.g. t0 and vdrift corrections of SDD
301 //
302 nrecElems = recArr->GetLast()+1;
303 if (nrecElems<2 || !recOpt.IsDigit()) {stopped = kTRUE; break;}
304 int idx = recOpt.Atoi();
305 UShort_t voluid = (idx<=kMaxITSSensID) ? GetModuleVolumeID(idx) : idx;
306 AliITSAlignMille2Module* mod = 0;
307 //
308 if (voluid>=kMinITSSupeModuleID) { // custom supermodule
309 for (int j=0; j<fNSuperModules; j++) {
310 if (voluid==GetSuperModule(j)->GetVolumeID()) {
311 mod = new AliITSAlignMille2Module(*GetSuperModule(j));
312 // the matrix might be updated in case some prealignment was applied, check
313 TGeoHMatrix* mup = AliGeomManager::GetMatrix(mod->GetName());
314 if (mup) *(mod->GetMatrix()) = *mup;
315 fMilleModule.AddAtAndExpand(mod,fNModules);
316 break;
317 }
318 }
319 }
320 else if (idx<=kMaxITSSensVID) {
321 mod = new AliITSAlignMille2Module(voluid);
322 fMilleModule.AddAtAndExpand(mod,fNModules);
323 }
324 if (!mod) {stopped = kTRUE; break;} // bad volid
325 //
326 // geometry variation settings
327 for (int i=0;i<AliITSAlignMille2Module::kMaxParGeom;i++) {
328 irec = i+2;
329 if (irec >= nrecElems) break;
330 recExt = recArr->At(irec)->GetName();
331 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
332 mod->SetFreeDOF(i, recExt.Atof() );
333 }
334 if (stopped) break;
335 //
336 // scaling factors for cluster errors
337 // first set default ones
338 for (int i=0;i<3;i++) mod->SetSigmaFactor(i, fSigmaFactor[i]);
339 for (int i=0;i<3;i++) {
340 irec = i+8;
341 if (irec >= nrecElems) break;
342 recExt = recArr->At(irec)->GetName();
343 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
344 mod->SetSigmaFactor(i, recExt.Atof() );
345 }
346 if (stopped) break;
347 //
348 mod->SetGeomParamsGlobal(fUseGlobalDelta);
349 // now comes special detectors treatment
350 if (mod->IsSDD()) {
351 double vl = 0;
352 if (nrecElems>11) {
353 recExt = recArr->At(11)->GetName();
354 if (recExt.IsFloat()) vl = recExt.Atof();
355 else {stopped = kTRUE; break;}
356 irec = 11;
357 }
358 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFT0,vl);
359 //
360 vl = 0;
361 if (nrecElems>12) {
362 recExt = recArr->At(12)->GetName();
363 if (recExt.IsFloat()) vl = recExt.Atof();
364 else {stopped = kTRUE; break;}
365 irec = 12;
366 }
367 mod->SetFreeDOF(AliITSAlignMille2Module::kDOFDV,vl);
368 }
369 //
370 mod->SetUniqueID(fNModules);
371 mod->EvaluateDOF();
372 fNModules++;
373 //
374 // now check if there are local constraints on this module
375 for (++irec;irec<nrecElems;irec++) {
376 recExt = recArr->At(irec)->GetName();
377 if (recExt.IsFloat()) {stopped=kTRUE;break;}
378 AliITSAlignMille2ConstrArray* cstr = (AliITSAlignMille2ConstrArray*)GetConstraint(recExt.Data());
379 if (!cstr) {
380 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
381 stopped=kTRUE;
382 break;
383 }
384 cstr->AddModule(mod);
385 }
386 if (stopped) break;
387 } // end while for loop over modules
388 if (stopped) break;
7b85e477 389 //
6526a72c 390 if (fNModules==0) {AliError("Failed to find any MODULE"); stopped = kTRUE; break;}
391 BuildHierarchy(); // preprocess loaded modules
7b85e477 392 //
6526a72c 393 // =========== 4: the rest may come in arbitrary order =======================================
394 rewind(pfc);
395 while ( (recArr=GetConfigRecord(pfc,recTitle="",recOpt,0))!=0 ) {
396 //
397 nrecElems = recArr->GetLast()+1;
398 //
399 // some simple flags -----------------------------------------------------------------------
400 //
401 if (recTitle == "SET_PSEUDO_PARENTS") SetAllowPseudoParents(kTRUE);
402 //
403 // some optional parameters ----------------------------------------------------------------
404 else if (recTitle == "SET_TRACK_FIT_METHOD") {
405 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
406 SetInitTrackParamsMeth(recOpt.Atoi());
7b85e477 407 }
6526a72c 408 //
409 else if (recTitle == "SET_MINPNT_TRA") {
410 if (recOpt.IsNull() || !recOpt.IsDigit() ) {stopped = kTRUE; break;}
411 fMinNPtsPerTrack = recOpt.Atoi();
7b85e477 412 }
7b85e477 413 //
6526a72c 414 else if (recTitle == "SET_NSTDDEV") {
415 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
416 fNStdDev = (Int_t)recOpt.Atof();
417 }
7b85e477 418 //
6526a72c 419 else if (recTitle == "SET_RESCUT_INIT") {
420 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
421 fResCutInitial = recOpt.Atof();
7b85e477 422 }
6526a72c 423 //
424 else if (recTitle == "SET_RESCUT_OTHER") {
425 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
426 fResCut = recOpt.Atof();
7b85e477 427 }
7b85e477 428 //
6526a72c 429 else if (recTitle == "SET_LOCALSIGMAFACTOR") { //-------------------------
430 for (irec=0;irec<3;irec++) if (nrecElems>irec+1) {
431 fSigmaFactor[irec] = ((TObjString*)recArr->At(irec+1))->GetString().Atof();
432 if (fSigmaFactor[irec]<=0.) stopped = kTRUE;
433 }
434 if (stopped) break;
435 }
7b85e477 436 //
6526a72c 437 else if (recTitle == "SET_STARTFAC") { //-------------------------
438 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
439 fStartFac = recOpt.Atof();
440 }
441 //
442 else if (recTitle == "SET_B_FIELD") { //-------------------------
443 if (recOpt.IsNull() || !recOpt.IsFloat() ) {stopped = kTRUE; break;}
444 fBField = recOpt.Atof();
445 if (fBField>0) {
446 fBOn = kTRUE;
447 fNLocal = 5; // helices
448 fRieman = new AliTrackFitterRieman();
449 }
450 else {
451 fBField = 0.0;
452 fBOn = kFALSE;
453 fNLocal = 4;
454 }
455 }
456 //
457 else if (recTitle == "SET_SPARSE_MATRIX") { // matrix solver type
458 //
459 AliMillePede2::SetGlobalMatSparse(kTRUE);
460 if (recOpt.IsNull()) continue;
461 // solver type and settings
462 if (recOpt == "MINRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolMinRes );
463 else if (recOpt == "FGMRES") AliMillePede2::SetIterSolverType( AliMinResSolve::kSolFGMRes );
464 else {stopped = kTRUE; break;}
465 //
466 if (nrecElems>=3) { // preconditioner type
467 recExt = recArr->At(2)->GetName();
468 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
469 AliMillePede2::SetMinResPrecondType( recExt.Atoi() );
470 }
471 //
472 if (nrecElems>=4) { // tolerance
473 recExt = recArr->At(3)->GetName();
474 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
475 AliMillePede2::SetMinResTol( recExt.Atof() );
476 }
477 //
478 if (nrecElems>=5) { // maxIter
479 recExt = recArr->At(4)->GetName();
480 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
481 AliMillePede2::SetMinResMaxIter( recExt.Atoi() );
482 }
483 }
484 //
485 else if (recTitle == "REQUIRE_POINT") { //-------------------------
486 // syntax: REQUIRE_POINT where ndet updw nreqpts
487 // where = LAYER or DETECTOR
488 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
489 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
490 // nreqpts = minimum number of points of that type
491 if (nrecElems>=5) {
492 recOpt.ToUpper();
493 int lr = ((TObjString*)recArr->At(2))->GetString().Atoi() - 1;
494 int hb = ((TObjString*)recArr->At(3))->GetString().Atoi();
495 int np = ((TObjString*)recArr->At(4))->GetString().Atoi();
496 fRequirePoints = kTRUE;
497 if (recOpt == "LAYER") {
498 if (lr<0 || lr>5) {stopped = kTRUE; break;}
499 if (hb>0) fNReqLayUp[lr] = np;
500 else if (hb<0) fNReqLayDown[lr] = np;
501 else fNReqLay[lr] = np;
502 }
503 else if (recOpt == "DETECTOR") {
504 if (lr<0 || lr>2) {stopped = kTRUE; break;}
505 if (hb>0) fNReqDetUp[lr] = np;
506 else if (hb<0) fNReqDetDown[lr] = np;
507 else fNReqDet[lr] = np;
508 }
509 else {stopped = kTRUE; break;}
510 }
511 else {stopped = kTRUE; break;}
512 }
513 //
514 // global constraints on the subunits/orphans
515 else if (recTitle == "CONSTRAINT_ORPHANS") { //------------------------
516 // expect CONSTRAINT_ORPHANS MEAN/MEDIAN Value parID0 ... parID1 ...
517 if (nrecElems<4) {stopped = kTRUE; break;}
518 recExt = recArr->At(2)->GetName();
519 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
520 double val = recExt.Atof();
521 UInt_t pattern = 0;
522 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
523 recExt = recArr->At(irec)->GetName();
524 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
525 pattern |= 0x1 << recExt.Atoi();
526 }
527 if (stopped) break;
528 if (recOpt == "MEAN") ConstrainOrphansMean(val,pattern);
529 else if (recOpt == "MEDIAN") ConstrainOrphansMedian(val,pattern);
530 else {stopped = kTRUE; break;}
531 }
532 //
533 else if (recTitle == "CONSTRAINT_SUBUNITS") { //------------------------
534 // expect ONSTRAINT_SUBUNITS MEAN/MEDIAN Value parID0 ... parID1 ... VolID1 ... VolIDn - VolIDm
535 if (nrecElems<5) {stopped = kTRUE; break;}
536 recExt = recArr->At(2)->GetName();
537 if (!recExt.IsFloat()) {stopped = kTRUE; break;}
538 double val = recExt.Atof();
539 UInt_t pattern = 0;
540 for (irec=3;irec<nrecElems;irec++) { // read params to constraint
541 recExt = recArr->At(irec)->GetName();
542 if (!recExt.IsDigit()) {stopped = kTRUE; break;}
543 int parid = recExt.Atoi();
544 if (parid<kMaxITSSensID) pattern |= 0x1 << recExt.Atoi();
545 else break; // list of params is over
546 }
547 if (stopped) break;
548 //
549 Bool_t meanC;
550 if (recOpt == "MEAN") meanC = kTRUE;
551 else if (recOpt == "MEDIAN") meanC = kFALSE;
552 else {stopped = kTRUE; break;}
553 //
554 int curID = -1;
555 int rangeStart = -1;
556 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
557 recExt = recArr->At(irec)->GetName();
558 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
559 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
560 else curID = recExt.Atoi();
561 //
562 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
563 // this was a range start or single
564 int start;
565 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
566 else start = curID; // create constraint either for single module (or 1st in the range)
567 for (int id=start;id<=curID;id++) {
568 int id0 = IsVIDDefined(id);
569 if (id0<0) {AliDebug(3,Form("Undefined module %d requested in the SubUnits constraint, skipping",id)); continue;}
570 if (meanC) ConstrainModuleSubUnitsMean(id0,val,pattern);
571 else ConstrainModuleSubUnitsMedian(id0,val,pattern);
572 }
573 }
574 if (rangeStart>=0) stopped = kTRUE; // unfinished range
575 if (stopped) break;
576 }
577 //
578 // association of modules with local constraints
579 else if (recTitle == "APPLY_CONSTRAINT") { //------------------------
580 // expect APPLY_CONSTRAINT NAME [NAME1...] [VolID1 ... VolIDn - VolIDm]
581 if (nrecElems<3) {stopped = kTRUE; break;}
582 int nmID0=-1,nmID1=-1;
583 for (irec=1;irec<nrecElems;irec++) { // find the range of constraint names
584 recExt = recArr->At(irec)->GetName();
585 if (recExt.IsFloat()) break;
586 // check if such a constraint was declared
587 if (!GetConstraint(recExt.Data())) {
588 AliInfo(Form("No Local constraint %s was declared",recExt.Data()));
589 stopped=kTRUE;
590 break;
591 }
592 if (nmID0<0) nmID0 = irec;
593 nmID1 = irec;
594 }
595 if (stopped) break;
596 //
597 if (irec>=nrecElems) {stopped = kTRUE; break;} // no modules provided
598 //
599 // now read the list of modules to constrain
600 int curID = -1;
601 int rangeStart = -1;
602 for (;irec<nrecElems;irec++) { // read modules to apply this constraint
603 recExt = recArr->At(irec)->GetName();
604 if (recExt == "-") {rangeStart = curID; continue;} // range is requested
605 else if (!recExt.IsDigit()) {stopped = kTRUE; break;}
606 else curID = recExt.Atoi();
607 //
608 if (curID<=kMaxITSSensID) curID = GetModuleVolumeID(curID);
609 //
610 // this was a range start or single
611 int start;
612 if (rangeStart>=0) {start = rangeStart+1; rangeStart=-1;} // continue the range
613 else start = curID; // create constraint either for single module (or 1st in the range)
614 for (int id=start;id<=curID;id++) {
615 AliITSAlignMille2Module *md = GetMilleModuleByVID(id);
616 if (!md) {AliDebug(3,Form("Undefined module %d requested in the Local constraint, skipping",id)); continue;}
617 for (int nmid=nmID0;nmid<=nmID1;nmid++)
618 ((AliITSAlignMille2ConstrArray*)GetConstraint(recArr->At(nmid)->GetName()))->AddModule(md);
619 }
620 }
621 if (rangeStart>=0) stopped = kTRUE; // unfinished range
622 if (stopped) break;
623 }
624 //
625 else continue; // already processed record
626 //
627 } // end of while loop 4 over the various params
7b85e477 628 //
6526a72c 629 break;
630 } // end of while(1) loop
7b85e477 631 //
632 fclose(pfc);
6526a72c 633 if (stopped) {
634 AliError(Form("Failed on record %s %s ...\n",recTitle.Data(),recOpt.Data()));
635 return -1;
636 }
637 //
638 fIsConfigured = kTRUE;
639 return 0;
640}
641
642//________________________________________________________________________________________________________
643void AliITSAlignMille2::BuildHierarchy()
644{
645 // build the hieararhy of the modules to align
7b85e477 646 //
6526a72c 647 if (!GetUseGlobalDelta() && PseudoParentsAllowed()) {
648 AliInfo("PseudoParents mode is allowed only when the deltas are global\n"
649 "Since Deltas are local, switching to NoPseudoParents");
650 SetAllowPseudoParents(kFALSE);
651 }
7b85e477 652 // set parent/child relationship for modules to align
6526a72c 653 AliInfo("Setting parent/child relationships\n");
7b85e477 654 //
6526a72c 655 // 1) child -> parent reference
656 for (int ipar=0;ipar<fNModules;ipar++) {
7b85e477 657 AliITSAlignMille2Module* parent = GetMilleModule(ipar);
6526a72c 658 if (parent->IsSensor()) continue; // sensor cannot be a parent
7b85e477 659 //
6526a72c 660 for (int icld=0;icld<fNModules;icld++) {
7b85e477 661 if (icld==ipar) continue;
662 AliITSAlignMille2Module* child = GetMilleModule(icld);
663 if (!child->BelongsTo(parent)) continue;
6526a72c 664 // child cannot have more sensors than the parent
665 if (child->GetNSensitiveVolumes() > parent->GetNSensitiveVolumes()) continue;
7b85e477 666 //
667 AliITSAlignMille2Module* parOld = child->GetParent();
6526a72c 668 // is this parent candidate closer than the old parent ?
7b85e477 669 if (parOld && parOld->GetNSensitiveVolumes()<parent->GetNSensitiveVolumes()) continue; // parOld is closer
670 child->SetParent(parent);
671 }
672 //
673 }
674 //
6526a72c 675 // add parent -> children reference
676 for (int icld=0;icld<fNModules;icld++) {
677 AliITSAlignMille2Module* child = GetMilleModule(icld);
678 AliITSAlignMille2Module* parent = child->GetParent();
679 if (parent) parent->AddChild(child);
680 }
681 //
7b85e477 682 // reorder the modules in such a way that parents come first
6526a72c 683 for (int icld=0;icld<fNModules;icld++) {
057304f3 684 AliITSAlignMille2Module* child = GetMilleModule(icld);
685 AliITSAlignMille2Module* parent;
6526a72c 686 while ( (parent=child->GetParent()) && (parent->GetUniqueID()>child->GetUniqueID()) ) {
7b85e477 687 // swap
688 fMilleModule[icld] = parent;
689 fMilleModule[parent->GetUniqueID()] = child;
690 child->SetUniqueID(parent->GetUniqueID());
691 parent->SetUniqueID(icld);
057304f3 692 child = parent;
7b85e477 693 }
057304f3 694 //
6526a72c 695 }
7b85e477 696 //
6526a72c 697 // Go over the child->parent chain and mark modules with explicitly provided sensors.
698 // If the sensors of the unit are explicitly declared, all undeclared sensors are
699 // suppresed in this unit.
700 for (int icld=fNModules;icld--;) {
7b85e477 701 AliITSAlignMille2Module* child = GetMilleModule(icld);
702 AliITSAlignMille2Module* parent = child->GetParent();
703 if (!parent) continue;
6526a72c 704 //
705 // check if this parent was already processed
706 if (!parent->AreSensorsProvided()) {
707 parent->DelSensitiveVolumes();
708 parent->SetSensorsProvided(kTRUE);
709 }
710 // reattach sensors to parent
711 for (int isc=child->GetNSensitiveVolumes();isc--;) {
712 UShort_t senVID = child->GetSensVolVolumeID(isc);
713 if (!parent->IsIn(senVID)) parent->AddSensitiveVolume(senVID);
7b85e477 714 }
715 }
716 //
7b85e477 717}
718
6526a72c 719// pepo
720//________________________________________________________________________________________________________
721void AliITSAlignMille2::SetCurrentModule(Int_t id)
722{
723 /// set the current supermodule
724 // new meaning
725 if (fMilleVersion>=2) {
726 fCurrentModule = GetMilleModule(id);
727 return;
728 }
729 // old meaning
730 if (fMilleVersion<=1) {
731 Int_t index=id;
732 /// set as current the SuperModule that contains the 'index' sens.vol.
733 if (index<0 || index>2197) {
734 AliInfo("index does not correspond to a sensitive volume!");
735 return;
736 }
737 UShort_t voluid=AliITSAlignMille2Module::GetVolumeIDFromIndex(index);
738 Int_t k=IsContained(voluid);
739 if (k>=0){
740 fCluster.SetVolumeID(voluid);
741 fCluster.SetXYZ(0,0,0);
742 InitModuleParams();
743 }
744 else
745 AliInfo(Form("module %d not defined\n",index));
746 }
747}
748// endpepo
749//________________________________________________________________________________________________________
7b85e477 750void AliITSAlignMille2::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts)
751{
752 // set minimum number of points in specific detector or layer
753 // where = LAYER or DETECTOR
754 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
755 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
756 // nreqpts = minimum number of points of that type
757 ndet--;
758 if (strstr(where,"LAYER")) {
759 if (ndet<0 || ndet>5) return;
760 if (updw>0) fNReqLayUp[ndet]=nreqpts;
761 else if (updw<0) fNReqLayDown[ndet]=nreqpts;
762 else fNReqLay[ndet]=nreqpts;
763 fRequirePoints=kTRUE;
764 }
765 else if (strstr(where,"DETECTOR")) {
766 if (ndet<0 || ndet>2) return;
767 if (updw>0) fNReqDetUp[ndet]=nreqpts;
768 else if (updw<0) fNReqDetDown[ndet]=nreqpts;
769 else fNReqDet[ndet]=nreqpts;
770 fRequirePoints=kTRUE;
771 }
772}
773
6526a72c 774//________________________________________________________________________________________________________
775Int_t AliITSAlignMille2::GetModuleIndex(const Char_t *symname)
776{
7b85e477 777 /// index from symname
778 if (!symname) return -1;
779 for (Int_t i=0;i<=kMaxITSSensID; i++) {
780 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
781 }
782 return -1;
783}
784
6526a72c 785//________________________________________________________________________________________________________
786Int_t AliITSAlignMille2::GetModuleIndex(UShort_t voluid)
787{
7b85e477 788 /// index from volume ID
789 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
790 if (lay<1|| lay>6) return -1;
791 Int_t idx=Int_t(voluid)-2048*lay;
792 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
793 for (Int_t ilay=1; ilay<lay; ilay++)
794 idx += AliGeomManager::LayerSize(ilay);
795 return idx;
796}
797
6526a72c 798//________________________________________________________________________________________________________
799UShort_t AliITSAlignMille2::GetModuleVolumeID(const Char_t *symname)
800{
7b85e477 801 /// volume ID from symname
802 /// works for sensitive volumes only
803 if (!symname) return 0;
804
805 for (UShort_t voluid=2000; voluid<13300; voluid++) {
806 Int_t modId;
807 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
808 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
809 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
810 }
811 }
812
813 return 0;
814}
815
6526a72c 816//________________________________________________________________________________________________________
817UShort_t AliITSAlignMille2::GetModuleVolumeID(Int_t index)
818{
7b85e477 819 /// volume ID from index
820 if (index<0) return 0;
821 if (index<2198)
822 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
823 else {
824 for (int i=0; i<fNSuperModules; i++) {
825 if (GetSuperModule(i)->GetIndex()==index) return GetSuperModule(i)->GetVolumeID();
826 }
827 }
828 return 0;
829}
830
6526a72c 831//________________________________________________________________________________________________________
832Int_t AliITSAlignMille2::InitGeometry()
833{
7b85e477 834 /// initialize geometry
6526a72c 835 AliInfo("Loading initial geometry");
7b85e477 836 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
837 fGeoManager = AliGeomManager::GetGeometry();
838 if (!fGeoManager) {
839 AliInfo("Couldn't initialize geometry");
6526a72c 840 return -1;
841 }
842 return 0;
843}
844
845//________________________________________________________________________________________________________
846Int_t AliITSAlignMille2::SetConstraintWrtRef(const char* reffname)
847{
848 // Load the global deltas from this file. The local gaussian constraints on some modules
849 // will be defined with respect to the deltas from this reference file, converted to local
850 // delta format. Note: conversion to local format requires reloading the geometry!
851 //
852 AliInfo(Form("Loading reference deltas for local constraints from %s",reffname));
853 if (!fGeoManager) return -1;
854 fConstrRefFileName = reffname;
855 if (fConstrRefFileName == "IDEAL") { // the reference is the ideal geometry, just create dummy reference array
856 fConstrRef = new TClonesArray("AliAlignObjParams",1);
857 return 0;
858 }
859 TFile *pref = TFile::Open(fConstrRefFileName.Data());
860 if (!pref->IsOpen()) return -2;
861 fConstrRef = (TClonesArray*)pref->Get("ITSAlignObjs");
862 pref->Close();
863 delete pref;
864 if (!fConstrRef) {
865 AliError(Form("Did not find reference prealignment deltas in %s",reffname));
866 return -1;
867 }
868 //
869 // we need ideal geometry to convert global deltas to local ones
870 if (fUsePreAlignment) {
871 AliError("The call of SetConstraintWrtRef must be done before application of the prealignment");
872 return -1;
873 }
874 //
875 AliInfo("Converting global reference deltas to local ones");
876 Int_t nprea = fConstrRef->GetEntriesFast();
877 for (int ix=0; ix<nprea; ix++) {
878 AliAlignObjParams *preo=(AliAlignObjParams*) fConstrRef->At(ix);
879 if (!preo->ApplyToGeometry()) return -1;
7b85e477 880 }
6526a72c 881 //
882 // now convert the global reference deltas to local ones
883 for (int i=fConstrRef->GetEntriesFast();i--;) {
884 AliAlignObjParams *preo = (AliAlignObjParams*)fConstrRef->At(i);
885 TGeoHMatrix * mupd = AliGeomManager::GetMatrix(preo->GetSymName());
886 if (!mupd) { // this is not alignable entry, need to look in the supermodules
887 for (int im=fNSuperModules;im--;) {
888 AliITSAlignMille2Module* mod = GetSuperModule(im);
889 if ( strcmp(mod->GetName(), preo->GetSymName()) ) continue;
890 mupd = mod->GetMatrix();
891 break;
892 }
893 if (!mupd) {
894 AliError(Form("Failed to find the volume for reference %s",preo->GetSymName()));
895 return -1;
896 }
897 }
898 TGeoHMatrix preMat;
899 preo->GetMatrix(preMat); // Delta_Glob
900 TGeoHMatrix tmpMat = *mupd; // Delta_Glob * Delta_Glob_Par * M
901 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
902 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
903 preo->SetMatrix(tmpMat); // local corrections
904 }
905 //
906 // we need to reload the geometry spoiled by this reference deltas...
907 delete fGeoManager;
908 AliInfo("Reloading initial geometry");
909 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
910 fGeoManager = AliGeomManager::GetGeometry();
911 return 0;
7b85e477 912}
913
6526a72c 914//________________________________________________________________________________________________________
915void AliITSAlignMille2::Init()
7b85e477 916{
6526a72c 917 //
918 if (fIsMilleInit) {
919 AliInfo("Millepede has been already initialized!");
920 return;
921 }
922 // range constraints in such a way that the childs are constrained before their parents
923 // orphan constraints come last
924 for (int ic=0;ic<GetNConstraints();ic++) {
925 for (int ic1=ic+1;ic1<GetNConstraints();ic1++) {
926 AliITSAlignMille2Constraint *cst0 = GetConstraint(ic);
927 AliITSAlignMille2Constraint *cst1 = GetConstraint(ic1);
928 if (cst0->GetModuleID()<cst1->GetModuleID()) {
929 // swap
930 fConstraints[ic] = cst1;
931 fConstraints[ic1] = cst0;
932 }
933 }
934 }
935 //
936 if (!GetUseGlobalDelta()) {
937 AliInfo("ATTENTION: The parameters are defined in the local frame, no check for degeneracy will be done");
938 for (int imd=fNModules;imd--;) {
939 AliITSAlignMille2Module* mod = GetMilleModule(imd);
940 int npar = mod->GetNParTot();
941 // the parameter may have max 1 free instance, otherwise the equations are underdefined
942 for (int ipar=0;ipar<npar;ipar++) {
943 if (!mod->IsFreeDOF(ipar)) continue;
944 mod->SetParOffset(ipar,fNGlobal++);
945 }
946 }
947 }
948 else {
949 // init millepede, decide which parameters are to be fitted explicitly
950 for (int imd=fNModules;imd--;) {
951 AliITSAlignMille2Module* mod = GetMilleModule(imd);
952 int npar = mod->GetNParTot();
953 // the parameter may have max 1 free instance, otherwise the equations are underdefined
954 for (int ipar=0;ipar<npar;ipar++) {
955 if (!mod->IsFreeDOF(ipar)) continue; // fixed
956 //
957 int nFreeInstances = 0;
958 //
959 AliITSAlignMille2Module* parent = mod;
960 Bool_t cstMeanMed=kFALSE,cstGauss=kFALSE;
961 //
962 Bool_t AddToFit = kFALSE;
963 // the parameter may be ommitted from explicit fit (if PseudoParentsAllowed is true) if
964 // 1) it is not explicitly constrained or its does not participate in Gaussian constraint
965 // 2) the same applies to all of its parents
966 // 3) it has at least 1 unconstrained direct child
967 while(parent) {
968 if (!parent->IsFreeDOF(ipar)) {parent = parent->GetParent(); continue;}
969 nFreeInstances++;
970 if (IsParModConstrained(parent,ipar, cstMeanMed, cstGauss)) nFreeInstances--;
971 if (cstGauss) AddToFit = kTRUE;
972 parent = parent->GetParent();
973 }
974 if (nFreeInstances>1) {
975 AliError(Form("Parameter#%d of module %s\nhas %d free instances in the "
976 "unconstrained parents\nSystem is undefined",ipar,mod->GetName(),nFreeInstances));
977 exit(1);
978 }
979 //
980 // i) Are PseudoParents allowed?
981 if (!PseudoParentsAllowed()) AddToFit = kTRUE;
982 // ii) check if this module has no child with such a free parameter. Since the order of this check
983 // goes from child to parent, by this moment such a parameter must have been already added
984 else if (!IsParModFamilyVaried(mod,ipar)) AddToFit = kTRUE; // no varied children at all
985 else if (!IsParFamilyFree(mod,ipar,1)) AddToFit = kTRUE; // no unconstrained direct children
986 // otherwise the value of this parameter can be extracted from simple contraint and the values of
987 // the relevant parameters of its children the fit is done. Hence it is not included
988 if (!AddToFit) continue;
989 //
990 // shall add this parameter to explicit fit
991 // printf("Adding %s %d -> %d\n",mod->GetName(), ipar, fNGlobal);
992 mod->SetParOffset(ipar,fNGlobal++);
993 }
994 }
995 }
996 //
997 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
998 fGlobalDerivatives = new Double_t[fNGlobal];
999 memset(fGlobalDerivatives,0,fNGlobal*sizeof(Double_t));
1000 //
1001 fMillepede->InitMille(fNGlobal,fNLocal,fNStdDev,fResCut,fResCutInitial);
7b85e477 1002 fIsMilleInit = kTRUE;
6526a72c 1003 //
1004 ResetLocalEquation();
1005 AliInfo("Parameters initialized to zero");
1006 //
7b85e477 1007 /// Fix non free parameters
1008 for (Int_t i=0; i<fNModules; i++) {
6526a72c 1009 AliITSAlignMille2Module* mod = GetMilleModule(i);
1010 for (Int_t j=0; j<mod->GetNParTot(); j++) {
1011 if (mod->GetParOffset(j)<0) continue; // not varied
1012 FixParameter(mod->GetParOffset(j),mod->GetParConstraint(j));
1013 fMillepede->SetParamGrID(i, mod->GetParOffset(j));
1014 }
7b85e477 1015 }
1016 //
1017 // Set iterations
6526a72c 1018 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
1019 //
7b85e477 1020}
1021
6526a72c 1022//________________________________________________________________________________________________________
1023void AliITSAlignMille2::AddConstraint(Double_t *par, Double_t value, Double_t sigma)
1024{
7b85e477 1025 /// Constrain equation defined by par to value
6526a72c 1026 if (!fIsMilleInit) Init();
1027 fMillepede->SetGlobalConstraint(par, value, sigma);
7b85e477 1028 AliInfo("Adding constraint");
1029}
1030
6526a72c 1031//________________________________________________________________________________________________________
1032void AliITSAlignMille2::InitGlobalParameters(Double_t *par)
1033{
7b85e477 1034 /// Initialize global parameters with par array
6526a72c 1035 if (!fIsMilleInit) Init();
7b85e477 1036 fMillepede->SetGlobalParameters(par);
1037 AliInfo("Init Global Parameters");
1038}
6526a72c 1039
1040//________________________________________________________________________________________________________
1041void AliITSAlignMille2::FixParameter(Int_t iPar, Double_t value)
1042{
7b85e477 1043 /// Parameter iPar is encourage to vary in [-value;value].
1044 /// If value == 0, parameter is fixed
1045 if (!fIsMilleInit) {
1046 AliInfo("Millepede has not been initialized!");
1047 return;
1048 }
1049 fMillepede->SetParSigma(iPar, value);
1050 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
1051}
1052
6526a72c 1053//________________________________________________________________________________________________________
7b85e477 1054void AliITSAlignMille2::ResetLocalEquation()
1055{
1056 /// Reset the derivative vectors
1057 for(int i=fNLocal;i--;) fLocalDerivatives[i] = 0.0;
6526a72c 1058 memset(fGlobalDerivatives, 0, fNGlobal*sizeof(double) );
7b85e477 1059}
1060
6526a72c 1061//________________________________________________________________________________________________________
7b85e477 1062Int_t AliITSAlignMille2::ApplyToGeometry()
1063{
1064 // apply starting realignment to ideal geometry
6526a72c 1065 AliInfo(Form("Using %s for prealignment",fPreAlignmentFileName.Data()));
7b85e477 1066 if (!fGeoManager) return -1;
6526a72c 1067 TFile *pref = TFile::Open(fPreAlignmentFileName.Data());
7b85e477 1068 if (!pref->IsOpen()) return -2;
1069 fPrealignment = (TClonesArray*)pref->Get("ITSAlignObjs");
1070 if (!fPrealignment) return -3;
1071 Int_t nprea = fPrealignment->GetEntriesFast();
1072 AliInfo(Form("Array of input misalignments with %d entries",nprea));
1073 //
1074 for (int ix=0; ix<nprea; ix++) {
1075 AliAlignObjParams *preo=(AliAlignObjParams*) fPrealignment->At(ix);
1076 Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(preo->GetVolUID());
1077 if (index>=0) {
1078 if (index>=fPreAlignQF.GetSize()) fPreAlignQF.Set(index+10);
1079 fPreAlignQF[index] = (int) preo->GetUniqueID()+1;
1080 }
1081 //TString nms = preo->GetSymName();
1082 //if (!nms.Contains("Ladder")) continue; //RRR
6526a72c 1083 //printf("Applying#%4d %s\n",ix,preo->GetSymName());
7b85e477 1084 if (!preo->ApplyToGeometry()) return -4;
1085 }
1086 //
1087 pref->Close();
1088 delete pref;
1089 //
1090 fUsePreAlignment = kTRUE;
1091 return 0;
1092}
1093
6526a72c 1094//________________________________________________________________________________________________________
7b85e477 1095Int_t AliITSAlignMille2::GetPreAlignmentQualityFactor(Int_t index) const
1096{
1097 if (!fUsePreAlignment || index<0 || index>=fPreAlignQF.GetSize()) return -1;
1098 return fPreAlignQF[index]-1;
1099}
1100
6526a72c 1101//________________________________________________________________________________________________________
1102AliTrackPointArray *AliITSAlignMille2::PrepareTrack(const AliTrackPointArray *atp)
1103{
7b85e477 1104 /// create a new AliTrackPointArray keeping only defined modules
1105 /// move points according to a given prealignment, if any
1106 /// sort alitrackpoints w.r.t. global Y direction, if selected
6526a72c 1107 const double kTiny = 1E-12;
1108 //
7b85e477 1109 AliTrackPointArray *atps=NULL;
1110 Int_t idx[20];
1111 Int_t npts=atp->GetNPoints();
1112
1113 /// checks if AliTrackPoints belong to defined modules
1114 Int_t ngoodpts=0;
1115 Int_t intidx[20];
1116
1117 for (int j=0; j<npts; j++) {
6526a72c 1118 intidx[j] = IsVIDContained(atp->GetVolumeID()[j]);
7b85e477 1119 if (intidx[j]>=0) ngoodpts++;
1120 }
6526a72c 1121 AliDebug(3,Form("Number of points in defined modules: %d out of %d",ngoodpts,npts));
7b85e477 1122
1123 // reject track if not enough points are left
1124 if (ngoodpts<fMinNPtsPerTrack) {
1125 AliInfo("Track with not enough points!");
1126 return NULL;
1127 }
1128 // >> RS
1129 AliTrackPoint p;
1130 // check points in specific places
1131 if (fRequirePoints) {
1132 Int_t nlayup[6],nlaydown[6],nlay[6];
1133 Int_t ndetup[3],ndetdown[3],ndet[3];
1134 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
1135 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
1136
1137 for (int i=0; i<npts; i++) {
1138 // skip not defined points
1139 if (intidx[i]<0) continue;
1140 Float_t xx=atp->GetX()[i];
1141 Float_t yy=atp->GetY()[i];
1142 Float_t r=TMath::Sqrt(xx*xx + yy*yy);
1143 int lay=-1;
1144 if (r<5) lay=0;
1145 else if (r>5 && r<10) lay=1;
1146 else if (r>10 && r<18) lay=2;
1147 else if (r>18 && r<30) lay=3;
1148 else if (r>30 && r<40) lay=4;
1149 else if (r>40) lay=5;
1150 if (lay<0) continue;
1151 int det=lay/2;
1152 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
1153
1154 if (yy>=0.0) { // UP point
1155 nlayup[lay]++;
1156 nlay[lay]++;
1157 ndetup[det]++;
1158 ndet[det]++;
1159 }
1160 else {
1161 nlaydown[lay]++;
1162 nlay[lay]++;
1163 ndetdown[det]++;
1164 ndet[det]++;
1165 }
1166 }
1167
1168 // checks minimum values
1169 Bool_t isok=kTRUE;
1170 for (Int_t j=0; j<6; j++) {
1171 if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE;
1172 if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE;
1173 if (nlay[j]<fNReqLay[j]) isok=kFALSE;
1174 }
1175 for (Int_t j=0; j<3; j++) {
1176 if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE;
1177 if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE;
1178 if (ndet[j]<fNReqDet[j]) isok=kFALSE;
1179 }
1180 if (!isok) {
1181 AliDebug(2,Form("Track does not meet all location point requirements!"));
1182 return NULL;
1183 }
1184 }
7b85e477 1185 // build a new track with (sorted) (prealigned) good points
6526a72c 1186 atps = (AliTrackPointArray*)fTrackBuff[ngoodpts-fMinNPtsPerTrack];
1187 if (!atps) {
1188 atps = new AliTrackPointArray(ngoodpts);
1189 fTrackBuff.AddAtAndExpand(atps,ngoodpts-fMinNPtsPerTrack);
1190 }
1191 //
1192 //
7b85e477 1193 for (int i=0; i<npts; i++) idx[i]=i;
1194 // sort track if required
6526a72c 1195 TMath::Sort(npts,atp->GetY(),idx); // sort descending...
1196 //
7b85e477 1197 Int_t npto=0;
1198 for (int i=0; i<npts; i++) {
1199 // skip not defined points
1200 if (intidx[idx[i]]<0) continue;
1201 atp->GetPoint(p,idx[i]);
1202
1203 // prealign point if required
1204 // get IDEAL matrix
6526a72c 1205 AliITSAlignMille2Module *mod = GetMilleModule(intidx[idx[i]]);
1206 TGeoHMatrix *svOrigMatrix = mod->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
7b85e477 1207 // get back real local coordinates: use OriginalGlobalMatrix because AliTrackPoints were written
1208 // with idel geometry
1209 Double_t pg[3],pl[3];
1210 pg[0]=p.GetX();
1211 pg[1]=p.GetY();
1212 pg[2]=p.GetZ();
1213 // printf("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1214 AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1215 svOrigMatrix->MasterToLocal(pg,pl);
1216
1217 AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]));
6526a72c 1218 //
1219 // this is a temporary code to extract the drift speed used for given point
1220 if (p.GetDriftTime()>0) { // RRR
1221 // calculate the drift speed
1222 int sid = AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());// - kSDDoffsID;
1223 fDriftTime0[npto] = fInitialRecSDD ? fInitialRecSDD->GetTimeZero(sid) : 0.;
1224 /*
1225 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(p.GetVolumeID());
1226 if (lay==3) fDriftTime0[npto] = pg[2]<0 ? 169.5 : 140.1;
1227 else if (lay==4) fDriftTime0[npto] = pg[2]<0 ? 158.3 : 139.0;
1228 else {
1229 AliError(Form("Strange layer %d for moduleID %d",lay,p.GetVolumeID()));
1230 exit(1);
1231 }
1232 */
1233 double tdif = p.GetDriftTime() - fDriftTime0[npto];
1234 if (tdif<=0) tdif = 1;
1235 double vdrift = (3.5085-TMath::Abs(pl[0]))/tdif;
1236 if (vdrift<0) vdrift = 0;
1237 //
1238 // TEMPORARY CORRECTION (if provided) -------------->>>
1239 if (fCorrectSDD) {
1240 float t0Upd = fCorrectSDD->GetTimeZero(sid);
1241 vdrift += fCorrectSDD->GetDeltaVDrift(sid);
1242 tdif = p.GetDriftTime() - t0Upd;
1243 // correct Xlocal
1244 pl[0] = TMath::Sign(3.5085 - vdrift*tdif,pl[0]);
1245 fDriftTime0[npto] = t0Upd;
1246 }
1247 // TEMPORARY CORRECTION (if provided) --------------<<<
1248 fDriftSpeed[npto] = TMath::Sign(vdrift,pl[0]);
1249 //
1250 /*
1251 printf("%d %+6.2f %+6.2f %+6.2f %+5.2f %+5.2f %+5.2f %+6.1f %+6.1f %+f %+f\n",
1252 p.GetVolumeID(),pg[0],pg[1],pg[2],pl[0],pl[1],pl[2],p.GetDriftTime(), fDriftTime0[npto], fDriftSpeed[npto],tdif);
1253 */
1254 }
7b85e477 1255
1256 // update covariance matrix
1257 TGeoHMatrix hcov;
1258 Double_t hcovel[9];
1259 hcovel[0]=double(p.GetCov()[0]);
1260 hcovel[1]=double(p.GetCov()[1]);
1261 hcovel[2]=double(p.GetCov()[2]);
1262 hcovel[3]=double(p.GetCov()[1]);
1263 hcovel[4]=double(p.GetCov()[3]);
1264 hcovel[5]=double(p.GetCov()[4]);
1265 hcovel[6]=double(p.GetCov()[2]);
1266 hcovel[7]=double(p.GetCov()[4]);
1267 hcovel[8]=double(p.GetCov()[5]);
1268 hcov.SetRotation(hcovel);
1269 // now rotate in local system
6526a72c 1270 // printf("\nErrMatGlob: before\n"); hcov.Print(""); //RRR
7b85e477 1271 hcov.Multiply(svOrigMatrix);
1272 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
1273 // now hcov is LOCAL COVARIANCE MATRIX
6526a72c 1274 // apply sigma scaling
1275 // printf("\nErrMatLoc: before\n"); hcov.Print(""); //RRR
1276 Double_t *hcovscl = hcov.GetRotationMatrix();
1277 // for (int ir=3;ir--;) for (int ic=3;ic--;) hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1278 // RS TEMPORARY: nullify non-diagonal elements and sigY
1279 hcovscl[5] = 0;
1280 for (int ir=3;ir--;) for (int ic=3;ic--;) {
1281 if (ir==ic) {
1282 if (TMath::Abs(hcovscl[ir*3+ic])<kTiny) hcovscl[ir*3+ic] = 0.;
1283 else hcovscl[ir*3+ic] *= mod->GetSigmaFactor(ir)*mod->GetSigmaFactor(ic); //RRR
1284 }
1285 else hcovscl[ir*3+ic] = 0;
1286 }
1287 //
1288 // printf("\nErrMatLoc: after\n"); hcov.Print(""); //RRR
1289 //
7b85e477 1290 if (fBug==1) {
1291 // correzione bug LAYER 5 SSD temporanea..
1292 int ssdidx=AliITSAlignMille2Module::GetIndexFromVolumeID(p.GetVolumeID());
1293 if (ssdidx>=500 && ssdidx<1248) {
1294 int ladder=(ssdidx-500)%22;
1295 if (ladder==18) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx+1));
1296 if (ladder==19) p.SetVolumeID(AliITSAlignMille2Module::GetVolumeIDFromIndex(ssdidx-1));
1297 }
1298 }
7b85e477 1299 /// get (evenctually prealigned) matrix of sens. vol.
6526a72c 1300 TGeoHMatrix *svMatrix = mod->GetSensitiveVolumeMatrix(p.GetVolumeID());
7b85e477 1301 // modify global coordinates according with pre-aligment
1302 svMatrix->LocalToMaster(pl,pg);
1303 // now rotate in local system
1304 hcov.Multiply(&svMatrix->Inverse());
1305 hcov.MultiplyLeft(svMatrix);
1306 // hcov is back in GLOBAL RF
6526a72c 1307 // cure once more
1308 for (int ir=3;ir--;) for (int ic=3;ic--;) if (TMath::Abs(hcovscl[ir*3+ic])<kTiny) hcovscl[ir*3+ic] = 0.;
1309 // printf("\nErrMatGlob: after\n"); hcov.Print(""); //RRR
1310 //
7b85e477 1311 Float_t pcov[6];
6526a72c 1312 pcov[0]=hcovscl[0];
1313 pcov[1]=hcovscl[1];
1314 pcov[2]=hcovscl[2];
1315 pcov[3]=hcovscl[4];
1316 pcov[4]=hcovscl[5];
1317 pcov[5]=hcovscl[8];
7b85e477 1318
1319 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1320 // printf("New Gl coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1321 AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
1322 atps->AddPoint(npto,&p);
6526a72c 1323 AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,atps->GetX()[npto],
1324 atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] ));
1325 // printf("Adding %d %d %f\n",npto, p.GetVolumeID(), p.GetY());
7b85e477 1326 npto++;
1327 }
1328
1329 return atps;
1330}
1331
6526a72c 1332//________________________________________________________________________________________________________
1333AliTrackPointArray *AliITSAlignMille2::SortTrack(const AliTrackPointArray *atp)
1334{
7b85e477 1335 /// sort alitrackpoints w.r.t. global Y direction
1336 AliTrackPointArray *atps=NULL;
1337 Int_t idx[20];
1338 Int_t npts=atp->GetNPoints();
1339 AliTrackPoint p;
1340 atps=new AliTrackPointArray(npts);
1341
1342 TMath::Sort(npts,atp->GetY(),idx);
1343
1344 for (int i=0; i<npts; i++) {
1345 atp->GetPoint(p,idx[i]);
1346 atps->AddPoint(i,&p);
1347 AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
1348 }
1349 return atps;
1350}
1351
6526a72c 1352//________________________________________________________________________________________________________
1353Int_t AliITSAlignMille2::GetCurrentLayer() const
1354{
1355 if (!fGeoManager) {
1356 AliInfo("ITS geometry not initialized!");
1357 return -1;
1358 }
1359 return (Int_t)AliGeomManager::VolUIDToLayer(fCluster.GetVolumeID());
1360}
7b85e477 1361
6526a72c 1362//________________________________________________________________________________________________________
1363Int_t AliITSAlignMille2::InitModuleParams()
1364{
7b85e477 1365 /// initialize geometry parameters for a given detector
1366 /// for current cluster (fCluster)
1367 /// fGlobalInitParam[] is set as:
1368 /// [tx,ty,tz,psi,theta,phi]
1369 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
1370 /// *** At the moment: using Raffalele's angles definition ***
1371 ///
1372 /// return 0 if success
6526a72c 1373 /// If module is found but has no parameters to vary, return 1
7b85e477 1374
1375 if (!fGeoManager) {
1376 AliInfo("ITS geometry not initialized!");
1377 return -1;
1378 }
1379
1380 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
1381
1382 // set the internal index (index in module list)
1383 UShort_t voluid=fCluster.GetVolumeID();
1384 //
1385 // IT IS VERY IMPORTANT: start from the end of the list, where the childs are located !!!
1386 Int_t k=fNModules-1;
1387 fCurrentModule = 0;
6526a72c 1388 // VERY IMPORTANT: if the sensors were explicitly provided, don't look in the supermodules
1389 while (k>=0 && ! (fCurrentModule=GetMilleModule(k))->IsIn(voluid)) k--;
7b85e477 1390 if (k<0) return -3;
6526a72c 1391 //
1392 /*
1393 // Check if the module has free params. If not, go over the parents
1394 AliITSAlignMille2Module* mdtmp = fCurrentModule;
1395 while (mdtmp && mdtmp->GetNParFree()==0) mdtmp = mdtmp->GetParent();
1396 if (!mdtmp) return 1; // nothing to vary here
1397 fCurrentModule = mdtmp;
1398 */
7b85e477 1399 //
1400 fModuleInitParam[0] = 0.0;
1401 fModuleInitParam[1] = 0.0;
1402 fModuleInitParam[2] = 0.0;
1403 fModuleInitParam[3] = 0.0; // psi (X)
1404 fModuleInitParam[4] = 0.0; // theta (Y)
1405 fModuleInitParam[5] = 0.0; // phi (Z)
6526a72c 1406 fModuleInitParam[6] = 0.0;
1407 fModuleInitParam[7] = 0.0;
7b85e477 1408 /// get (evenctually prealigned) matrix of sens. vol.
1409 TGeoHMatrix *svMatrix = fCurrentModule->GetSensitiveVolumeMatrix(voluid);
1410
1411 fMeasGlo[0] = fCluster.GetX();
1412 fMeasGlo[1] = fCluster.GetY();
1413 fMeasGlo[2] = fCluster.GetZ();
1414 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
1415 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
1416
1417 // set stdev from cluster
1418 TGeoHMatrix hcov;
1419 Double_t hcovel[9];
1420 hcovel[0]=double(fCluster.GetCov()[0]);
1421 hcovel[1]=double(fCluster.GetCov()[1]);
1422 hcovel[2]=double(fCluster.GetCov()[2]);
1423 hcovel[3]=double(fCluster.GetCov()[1]);
1424 hcovel[4]=double(fCluster.GetCov()[3]);
1425 hcovel[5]=double(fCluster.GetCov()[4]);
1426 hcovel[6]=double(fCluster.GetCov()[2]);
1427 hcovel[7]=double(fCluster.GetCov()[4]);
1428 hcovel[8]=double(fCluster.GetCov()[5]);
1429 hcov.SetRotation(hcovel);
1430 // now rotate in local system
1431 hcov.Multiply(svMatrix);
1432 hcov.MultiplyLeft(&svMatrix->Inverse());
6526a72c 1433 //
7b85e477 1434 // set local sigmas
1435 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
1436 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4])); // RS
1437 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
1438
1439 // set minimum value for SigmaLoc to 10 micron
1440 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
1441 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
6526a72c 1442 //
1443 /* RRR the rescaling is moved to PrepareTrack
7b85e477 1444 // multiply local sigmas by global and module specific factor
1445 for (int i=3;i--;) fSigmaLoc[i] *= fSigmaFactor[i]*fCurrentModule->GetSigmaFactor(i);
1446 //
6526a72c 1447 */
7b85e477 1448 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
1449
1450 return 0;
1451}
1452
6526a72c 1453//________________________________________________________________________________________________________
7b85e477 1454void AliITSAlignMille2::Print(Option_t*) const
1455{
1456 ///
1457 printf("*** AliMillepede for ITS ***\n");
1458 printf(" Number of defined super modules: %d\n",fNModules);
1459 printf(" Obtained parameters refer to %s Deltas\n",fUseGlobalDelta ? "GLOBAL":"LOCAL");
1460 //
1461 if (fGeoManager)
1462 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
1463 else
1464 printf(" geometry not loaded\n");
1465 //
1466 if (fUsePreAlignment)
1467 printf(" using prealignment from %s \n",fPreAlignmentFileName.Data());
1468 else
1469 printf(" prealignment not used\n");
1470 //
1471 //
1472 if (fBOn)
1473 printf(" B Field set to %f T - using Riemann's helices\n",fBField);
1474 else
1475 printf(" B Field OFF - using straight lines \n");
1476 //
1477 if (fRequirePoints) printf(" Required points in tracks:\n");
1478 for (Int_t i=0; i<6; i++) {
1479 if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
1480 if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
1481 if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]);
1482 }
1483 for (Int_t i=0; i<3; i++) {
1484 if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
1485 if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
1486 if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]);
1487 }
1488 //
1489 printf("\n Millepede configuration parameters:\n");
7b85e477 1490 printf(" init value for chi2 cut : %.4f\n",fStartFac);
1491 printf(" first iteration cut value : %.4f\n",fResCutInitial);
1492 printf(" other iterations cut value : %.4f\n",fResCut);
1493 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
6526a72c 1494 printf(" def.scaling for local sigmas : %.4f %.4f %.4f\n",fSigmaFactor[0],fSigmaFactor[1],fSigmaFactor[2]);
7b85e477 1495
1496 printf("List of defined modules:\n");
1497 printf(" intidx\tindex\tvoluid\tname\n");
1498 for (int i=0; i<fNModules; i++) {
1499 AliITSAlignMille2Module* md = GetMilleModule(i);
1500 printf(" %d\t%d\t%d\t%s\n",i,md->GetIndex(),md->GetVolumeID(),md->GetName());
1501 }
1502}
1503
6526a72c 1504//________________________________________________________________________________________________________
7b85e477 1505AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleByVID(UShort_t voluid) const
1506{
6526a72c 1507 // return pointer to a defined supermodule
7b85e477 1508 // return NULL if error
6526a72c 1509 Int_t i=IsVIDDefined(voluid);
7b85e477 1510 if (i<0) return NULL;
1511 return GetMilleModule(i);
1512}
1513
6526a72c 1514//________________________________________________________________________________________________________
1515AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleBySymName(const Char_t* symname) const
7b85e477 1516{
6526a72c 1517 // return pointer to a defined supermodule
1518 // return NULL if error
1519 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1520 if (vid>0) return GetMilleModuleByVID(vid);
1521 else { // this is not alignable module, need to look within defined supermodules
1522 int i = IsSymDefined(symname);
1523 if (i>=0) return GetMilleModule(i);
1524 }
1525 return 0;
1526}
1527
1528//________________________________________________________________________________________________________
1529AliITSAlignMille2Module *AliITSAlignMille2::GetMilleModuleIfContained(const Char_t* symname) const
1530{
1531 // return pointer to a defined/contained supermodule
1532 // return NULL otherwise
1533 int i = IsSymContained(symname);
1534 return i<0 ? 0 : GetMilleModule(i);
1535}
1536
1537//________________________________________________________________________________________________________
1538AliAlignObjParams* AliITSAlignMille2::GetPrealignedObject(const Char_t* symname) const
1539{
1540 if (!fPrealignment) return 0;
1541 for (int ipre=fPrealignment->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1542 AliAlignObjParams* preob = (AliAlignObjParams*)fPrealignment->At(ipre);
1543 if (!strcmp(preob->GetSymName(),symname)) return preob;
1544 }
1545 return 0;
1546}
1547
1548//________________________________________________________________________________________________________
1549AliAlignObjParams* AliITSAlignMille2::GetConstrRefObject(const Char_t* symname) const
1550{
1551 if (!fConstrRef) return 0;
1552 for (int ipre=fConstrRef->GetEntriesFast();ipre--;) { // was the corresponding object prealigned?
1553 AliAlignObjParams* preob = (AliAlignObjParams*)fConstrRef->At(ipre);
1554 if (!strcmp(preob->GetSymName(),symname)) return preob;
1555 }
1556 return 0;
1557}
1558
1559//________________________________________________________________________________________________________
1560Bool_t AliITSAlignMille2::InitRiemanFit()
1561{
1562 // Initialize Riemann Fitter for current track
1563 // return kFALSE if error
7b85e477 1564
1565 if (!fBOn) return kFALSE;
1566
1567 Int_t npts=0;
1568 AliTrackPoint ap;
1569 npts = fTrack->GetNPoints();
1570 AliDebug(3,Form("Fitting track with %d points",npts));
1571
1572 fRieman->Reset();
1573 fRieman->SetTrackPointArray(fTrack);
1574
1575 TArrayI ai(npts);
1576 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
1577
1578 // fit track with 5 params in his own tracking-rotated reference system
1579 // xc = -p[1]/p[0];
1580 // yc = 1/p[0];
1581 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
1582 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
1583 return kFALSE;
1584 }
1585
1586 for (int i=0; i<5; i++)
1587 fLocalInitParam[i] = fRieman->GetParam()[i];
1588
1589 return kTRUE;
1590}
1591
6526a72c 1592//________________________________________________________________________________________________________
7b85e477 1593Bool_t fullErr2D = kTRUE;
1594
1595void trackFit2D(Int_t &, Double_t *, double &chi2, double *par, int)
1596{
1597 const double kTiny = 1.e-14;
1598 chi2 = 0;
1599 static AliTrackPoint pnt;
1600 //
1601 enum {kAX,kAZ,kBX,kBZ};
1602 enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5};
1603 //
1604 AliITSAlignMille2* alig = AliITSAlignMille2::GetInstance();
1605 AliTrackPointArray* track = alig->GetCurrentTrack();
1606 //
1607 int npts = track->GetNPoints();
1608 for (int ip=0;ip<npts;ip++) {
1609 track->GetPoint(pnt,ip);
1610 const float *cov = pnt.GetCov();
1611 double y = pnt.GetY();
1612 double dx = pnt.GetX() - (par[kAX]+y*par[kBX]);
1613 double dz = pnt.GetZ() - (par[kAZ]+y*par[kBZ]);
1614 double xxe = cov[kXX];
1615 double zze = cov[kZZ];
1616 double xze = cov[kXZ];
1617 //
1618 if (fullErr2D) {
1619 xxe += par[kBX]*par[kBX]*cov[kYY]-2.*par[kBX]*cov[kXY];
1620 zze += par[kBZ]*par[kBZ]*cov[kYY]-2.*par[kBZ]*cov[kZY];
1621 xze += par[kBX]*par[kBZ]*cov[kYY]-cov[kYZ]*par[kBZ]-cov[kXY]*par[kBX];
1622 }
1623 //
1624 double det = xxe*zze - xze*xze;
1625 if (det<kTiny) {
1626 printf("Negative diag. error (det=%+e) |sxx:%+e szz:%+e sxz:%+e| bx:%+e bz:%+e|\n"
1627 "Discarding correlation term\n",det,xxe,zze,xze,par[kBX],par[kBZ]);
1628 xxe = cov[kXX];
1629 zze = cov[kZZ];
1630 xze = cov[kXZ];
1631 fullErr2D = kFALSE;
1632 }
1633 double xxeI = zze/det;
1634 double zzeI = xxe/det;
1635 double xzeI =-xze/det;
1636 //
1637 chi2 += dx*dx*xxeI + dz*dz*zzeI + 2.*dx*dz*xzeI;
1638 //
1639 // printf("%d | %+e %+e %+e %+e %+e -> %+e\n",ip,dx,dz,xxeI,zzeI,xzeI, chi2);
1640 }
1641 //
1642}
1643
6526a72c 1644//________________________________________________________________________________________________________
7b85e477 1645void AliITSAlignMille2::InitTrackParams(int meth)
1646{
1647 /// initialize local parameters with different methods
1648 /// for current track (fTrack)
1649 Int_t npts=0;
1650 AliTrackPoint ap;
1651 double sX=0,sXY=0,sZ=0,sZY=0,sY=0,sYY=0,det=0;
1652 // simple linear interpolation
1653 // get local starting parameters (to be substituted by ESD track parms)
1654 // local parms (fLocalInitParam[]) are:
1655 // [0] = global x coord. of straight line intersection at y=0 plane
1656 // [1] = global z coord. of straight line intersection at y=0 plane
1657 // [2] = px/py
1658 // [3] = pz/py
1659 // test #1: linear fit in x(y) and z(y)
1660 npts = fTrack->GetNPoints();
1661 AliDebug(3,Form("*** initializing track with %d points ***",npts));
1662 for (int i=npts;i--;) {
1663 sY += fTrack->GetY()[i];
1664 sYY += fTrack->GetY()[i]*fTrack->GetY()[i];
1665 sX += fTrack->GetX()[i];
1666 sXY += fTrack->GetX()[i]*fTrack->GetY()[i];
1667 sZ += fTrack->GetZ()[i];
1668 sZY += fTrack->GetZ()[i]*fTrack->GetY()[i];
1669 }
1670 det = sYY*npts-sY*sY;
1671 if (det==0) det = 1E-20;
1672 fLocalInitParam[0] = (sX*sYY-sY*sXY)/det;
1673 fLocalInitParam[2] = (sXY*npts-sY*sX)/det;
1674 //
1675 fLocalInitParam[1] = (sZ*sYY-sY*sZY)/det;
1676 fLocalInitParam[3] = (sZY*npts-sY*sZ)/det;
7b85e477 1677 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f ugx = %f\n",fLocalInitParam[0],fLocalInitParam[2]));
1678 //
1679 if (meth==1) return;
1680 //
1681 // perform full fit accounting for cov.matrix
1682 static TVirtualFitter *minuit = 0;
1683 static Double_t step[5] = {1E-3,1E-3,1E-4,1E-4,1E-5};
1684 static Double_t arglist[10];
1685 //
1686 if (!minuit) {
1687 minuit = TVirtualFitter::Fitter(0,4);
1688 minuit->SetFCN(trackFit2D);
1689 arglist[0] = 1;
1690 minuit->ExecuteCommand("SET ERR",arglist, 1);
1691 //
1692 arglist[0] = -1;
1693 minuit->ExecuteCommand("SET PRINT",arglist,1);
1694 //
1695 }
1696 //
1697 minuit->SetParameter(0, "ax", fLocalInitParam[0], step[0], 0,0);
1698 minuit->SetParameter(1, "az", fLocalInitParam[1], step[1], 0,0);
1699 minuit->SetParameter(2, "bx", fLocalInitParam[2], step[2], 0,0);
1700 minuit->SetParameter(3, "bz", fLocalInitParam[3], step[3], 0,0);
1701 //
1702 arglist[0] = 1000; // number of function calls
1703 arglist[1] = 0.001; // tolerance
6526a72c 1704 fullErr2D = kFALSE;//kTRUE;
7b85e477 1705 minuit->ExecuteCommand("MIGRAD",arglist,2);
1706 //
1707 for (int i=0;i<4;i++) fLocalInitParam[i] = minuit->GetParameter(i);
1708 for (int i=0;i<4;i++) for (int j=0;j<4;j++) fLocalInitParEr[i][j] = minuit->GetCovarianceMatrixElement(i,j);
1709 //
1710}
1711
6526a72c 1712//________________________________________________________________________________________________________
1713Int_t AliITSAlignMille2::IsSymDefined(const Char_t* symname) const
1714{
1715 // checks if supermodule with this symname is defined and return the internal index
1716 // return -1 if not.
1717 for (int k=fNModules;k--;) if (!strcmp(symname,GetMilleModule(k)->GetName())) return k;
1718 return -1;
1719}
1720
1721//________________________________________________________________________________________________________
1722Int_t AliITSAlignMille2::IsSymContained(const Char_t* symname) const
1723{
1724 // checks if module with this symname is defined and return the internal index
1725 // return -1 if not.
1726 UShort_t vid = AliITSAlignMille2Module::GetVolumeIDFromSymname(symname);
1727 if (vid>0) return IsVIDContained(vid);
1728 // only sensors have real vid, but maybe we have a supermodule with fake vid?
1729 // IMPORTANT: always start from the end to start from the sensors
1730 return IsSymDefined(symname);
1731}
7b85e477 1732
6526a72c 1733//________________________________________________________________________________________________________
1734Int_t AliITSAlignMille2::IsVIDDefined(UShort_t voluid) const
7b85e477 1735{
1736 // checks if supermodule 'voluid' is defined and return the internal index
6526a72c 1737 // return -1 if not.
1738 for (int k=fNModules;k--;) if (voluid==GetMilleModule(k)->GetVolumeID()) return k;
1739 return -1;
7b85e477 1740}
1741
6526a72c 1742//________________________________________________________________________________________________________
1743Int_t AliITSAlignMille2::IsVIDContained(UShort_t voluid) const
7b85e477 1744{
6526a72c 1745 // checks if the sensitive module 'voluid' is contained inside a supermodule
1746 // and return the internal index of the last identified supermodule
7b85e477 1747 // return -1 if error
6526a72c 1748 // IMPORTANT: always start from the end to start from the sensors
7b85e477 1749 if (AliITSAlignMille2Module::GetIndexFromVolumeID(voluid)<0) return -1;
6526a72c 1750 for (int k=fNModules;k--;) if (GetMilleModule(k)->IsIn(voluid)) return k;
1751 return -1;
7b85e477 1752}
1753
6526a72c 1754//________________________________________________________________________________________________________
1755Int_t AliITSAlignMille2::CheckCurrentTrack()
7b85e477 1756{
7b85e477 1757 /// checks if AliTrackPoints belongs to defined modules
1758 /// return number of good poins
1759 /// return 0 if not enough points
1760
1761 Int_t npts = fTrack->GetNPoints();
1762 Int_t ngoodpts=0;
1763 // debug points
6526a72c 1764 for (int j=0; j<npts; j++) if (IsVIDContained(fTrack->GetVolumeID()[j])>=0) ngoodpts++;
1765 //
7b85e477 1766 if (ngoodpts<fMinNPtsPerTrack) return 0;
1767
1768 return ngoodpts;
1769}
1770
6526a72c 1771//________________________________________________________________________________________________________
1772Int_t AliITSAlignMille2::ProcessTrack(const AliTrackPointArray *track)
1773{
7b85e477 1774 /// Process track; Loop over hits and set local equations
1775 /// here 'track' is a AliTrackPointArray
1776 /// return 0 if success;
1777
6526a72c 1778 if (!fIsMilleInit) Init();
1779 //
7b85e477 1780 Int_t npts = track->GetNPoints();
1781 AliDebug(2,Form("*** Input track with %d points ***",npts));
1782
1783 // preprocessing of the input track: keep only points in defined volumes,
1784 // move points if prealignment is set, sort by Yglo if required
1785
1786 fTrack=PrepareTrack(track);
1787 if (!fTrack) return -1;
1788
1789 npts = fTrack->GetNPoints();
6526a72c 1790 if (npts>kMaxPoints) {
1791 AliError(Form("Compiled with kMaxPoints=%d, current track has %d points",kMaxPoints,npts));
1792 }
7b85e477 1793 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
1794
1795 if (!fBOn) { // straight lines
1796 // set local starting parameters (to be substituted by ESD track parms)
1797 // local parms (fLocalInitParam[]) are:
1798 // [0] = global x coord. of straight line intersection at y=0 plane
1799 // [1] = global z coord. of straight line intersection at y=0 plane
1800 // [2] = px/py
1801 // [3] = pz/py
1802 InitTrackParams(fInitTrackParamsMeth);
1803 }
1804 else {
1805 // local parms (fLocalInitParam[]) are the Riemann Fitter params
1806 if (!InitRiemanFit()) {
1807 AliInfo("Riemann fit failed! skipping this track...");
7b85e477 1808 fTrack=NULL;
1809 return -5;
1810 }
1811 }
1812
1813 Int_t nloceq=0;
6526a72c 1814 Int_t ngloeq=0;
1815 static Mille2Data md[kMaxPoints];
1816 //
7b85e477 1817 for (Int_t ipt=0; ipt<npts; ipt++) {
1818 fTrack->GetPoint(fCluster,ipt);
6526a72c 1819 fCluster.SetUniqueID(ipt);
7b85e477 1820 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
1821
1822 // set geometry parameters for the the current module
1823 if (InitModuleParams()) continue;
1824 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n",
1825 track->GetVolumeID()[ipt], fCurrentModule->GetIndex(),
1826 fCurrentModule->GetUniqueID(), AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
1827 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
6526a72c 1828 int res = AddLocalEquation(md[nloceq]);
1829 if (res<0) {fTotBadLocEqPoints++; nloceq = 0; break;}
1830 else if (res==0) nloceq++;
1831 else {nloceq++; ngloeq++;}
7b85e477 1832 } // end loop over points
1833 //
7b85e477 1834 fTrack=NULL;
6526a72c 1835 // not enough good points?
1836 if (nloceq<fMinNPtsPerTrack || ngloeq<1) return -1;
7b85e477 1837 //
1838 // finally send local equations to millepede
1839 SetLocalEquations(md,nloceq);
1840 fMillepede->SaveRecordData(); // RRR
1841 //
1842 return 0;
1843}
1844
6526a72c 1845//________________________________________________________________________________________________________
1846Int_t AliITSAlignMille2::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar)
1847{
7b85e477 1848 /// calculate track intersection point in local coordinates
1849 /// according with a given set of parameters (local(4) and global(6))
1850 /// and fill fPintLoc/Glo
1851 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
1852 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
1853 /// return 0 if success
1854
1855 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]));
1856 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]));
1857
1858
1859 // prepare the TGeoHMatrix
6526a72c 1860 TGeoHMatrix *tempHMat = fCurrentModule->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar,
1861 !fUseGlobalDelta);
1862 if (!tempHMat) return -1;
7b85e477 1863
1864 Double_t v0g[3]; // vector with straight line direction in global coord.
1865 Double_t p0g[3]; // point of the straight line (glo)
1866
1867 if (fBOn) { // B FIELD!
1868 AliTrackPoint prf;
1869 for (int ip=0; ip<5; ip++)
1870 fRieman->SetParam(ip,lpar[ip]);
1871
1872 if (!fRieman->GetPCA(fCluster,prf)) {
1873 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
1874 return -3;
1875 }
1876 // now determine straight line passing tangent to fit curve at prf
1877 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
1878 // mo' P1=(X,Y,Z)_glo_prf
1879 // => (x,y,Z)_trk_prf ruotando di alpha...
1880 Double_t alpha=fRieman->GetAlpha();
1881 Double_t x1g = prf.GetX();
1882 Double_t y1g = prf.GetY();
1883 Double_t z1g = prf.GetZ();
1884 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
1885 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
1886 Double_t z1t = z1g;
1887
1888 Double_t x2t = x1t+1.0;
1889 Double_t y2t = y1t+fRieman->GetDYat(x1t);
1890 Double_t z2t = z1t+fRieman->GetDZat(x1t);
1891 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
1892 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
1893 Double_t z2g = z2t;
1894
1895 AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi()));
1896 AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
1897 AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
1898
1899 if (TMath::Abs(y2g-y1g)<1e-15) {
1900 AliInfo("DeltaY=0! Cannot proceed...");
1901 return -1;
1902 }
1903 // ugx,1,ugz
1904 v0g[0] = (x2g-x1g)/(y2g-y1g);
1905 v0g[1]=1.0;
1906 v0g[2] = (z2g-z1g)/(y2g-y1g);
1907
1908 // point: just keep prf
1909 p0g[0]=x1g;
1910 p0g[1]=y1g;
1911 p0g[2]=z1g;
1912 }
1913 else { // staight line
1914 // vector of initial straight line direction in glob. coord
1915 v0g[0]=lpar[2];
1916 v0g[1]=1.0;
1917 v0g[2]=lpar[3];
1918
1919 // intercept in yg=0 plane in glob coord
1920 p0g[0]=lpar[0];
1921 p0g[1]=0.0;
1922 p0g[2]=lpar[1];
1923 }
1924 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]));
1925
1926 // same in local coord.
1927 Double_t p0l[3],v0l[3];
6526a72c 1928 tempHMat->MasterToLocalVect(v0g,v0l);
1929 tempHMat->MasterToLocal(p0g,p0l);
7b85e477 1930
1931 if (TMath::Abs(v0l[1])<1e-15) {
1932 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
1933 return -1;
1934 }
1935
1936 // local intersection point
1937 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
1938 fPintLoc[1] = 0;
1939 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
1940
1941 // global intersection point
6526a72c 1942 tempHMat->LocalToMaster(fPintLoc,fPintGlo);
7b85e477 1943 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]));
1944
1945 return 0;
1946}
1947
6526a72c 1948//________________________________________________________________________________________________________
1949Int_t AliITSAlignMille2::CalcDerivatives(Int_t paridx, Bool_t islpar)
1950{
7b85e477 1951 /// calculate numerically (ROOT's style) the derivatives for
1952 /// local X intersection and local Z intersection
1953 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
1954 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
1955 /// return 0 if success
1956
1957 // copy initial parameters
6526a72c 1958 Double_t lpar[kNLocal];
1959 Double_t gpar[kNParCh];
1960 Double_t *derivative;
1961 for (Int_t i=0; i<kNLocal; i++) lpar[i]=fLocalInitParam[i];
1962 for (Int_t i=0; i<kNParCh; i++) gpar[i]=fModuleInitParam[i];
7b85e477 1963
1964 // trial with fixed dpar...
6526a72c 1965 Double_t dpar = 0.;
7b85e477 1966
1967 if (islpar) { // track parameters
1968 //dpar=fLocalInitParam[paridx]*0.001;
1969 // set minimum dpar
6526a72c 1970 derivative = fDerivativeLoc[paridx];
7b85e477 1971 if (!fBOn) {
6526a72c 1972 if (paridx<3) dpar=1.0e-4; // translations
7b85e477 1973 else dpar=1.0e-6; // direction
1974 }
1975 else { // B Field
1976 // pepo: proviamo con 1/1000, poi evenctually 1/100...
1977 Double_t dfrac=0.01;
1978 switch(paridx) {
1979 case 0:
1980 // RMS cosmics: 1e-4
1981 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1982 break;
1983 case 1:
1984 // RMS cosmics: 0.2
1985 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1986 break;
1987 case 2:
1988 // RMS cosmics: 9
1989 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1990 break;
1991 case 3:
1992 // RMS cosmics: 7
1993 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1994 break;
1995 case 4:
1996 // RMS cosmics: 0.3
1997 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1998 break;
1999 }
2000 }
2001 }
2002 else { // alignment global parameters
6526a72c 2003 derivative = fDerivativeGlo[paridx];
7b85e477 2004 //dpar=fModuleInitParam[paridx]*0.001;
2005 if (paridx<3) dpar=1.0e-4; // translations
2006 else dpar=1.0e-2; // angles
2007 }
2008
2009 AliDebug(3,Form("+++ using dpar=%g",dpar));
2010
2011 // calculate derivative ROOT's like:
2012 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
2013 Double_t pintl1[3]; // f(x-h)
2014 Double_t pintl2[3]; // f(x-h/2)
2015 Double_t pintl3[3]; // f(x+h/2)
2016 Double_t pintl4[3]; // f(x+h)
2017
2018 // first values
2019 if (islpar) lpar[paridx] -= dpar;
2020 else gpar[paridx] -= dpar;
2021 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2022 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
2023
2024 // second values
2025 if (islpar) lpar[paridx] += dpar/2;
2026 else gpar[paridx] += dpar/2;
2027 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2028 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
2029
2030 // third values
2031 if (islpar) lpar[paridx] += dpar;
2032 else gpar[paridx] += dpar;
2033 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2034 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
2035
2036 // fourth values
2037 if (islpar) lpar[paridx] += dpar/2;
2038 else gpar[paridx] += dpar/2;
2039 if (CalcIntersectionPoint(lpar, gpar)) return -2;
2040 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
2041
2042 Double_t h2 = 1./(2.*dpar);
2043 Double_t d0 = pintl4[0]-pintl1[0];
2044 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
6526a72c 2045 derivative[0] = h2*(4*d2 - d0)/3.;
2046 if (TMath::Abs(derivative[0]) < 1.0e-9) derivative[0] = 0.0;
7b85e477 2047
2048 d0 = pintl4[2]-pintl1[2];
2049 d2 = 2.*(pintl3[2]-pintl2[2]);
6526a72c 2050 derivative[2] = h2*(4*d2 - d0)/3.;
2051 if (TMath::Abs(derivative[2]) < 1.0e-9) derivative[2]=0.0;
7b85e477 2052
2053 AliDebug(3,Form("\n+++ derivatives +++ \n"));
6526a72c 2054 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",derivative[0]));
2055 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",derivative[2]));
7b85e477 2056
2057 return 0;
2058}
2059
6526a72c 2060//________________________________________________________________________________________________________
2061Int_t AliITSAlignMille2::AddLocalEquation(Mille2Data &m)
2062{
7b85e477 2063 /// Define local equation for current cluster in X and Z coor.
2064 /// and store them to memory
6526a72c 2065 /// return -1 in case of failure to build some equation
2066 /// 0 if no free global parameters were found but local eq is built
2067 /// 1 if both local and global eqs are built
2068 //
7b85e477 2069 // store first interaction point
6526a72c 2070 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -1;
7b85e477 2071 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
2072 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
2073
2074 // calculate local derivatives numerically
2075 Bool_t zeroX = kTRUE;
2076 Bool_t zeroZ = kTRUE;
2077 //
2078 for (Int_t i=0; i<fNLocal; i++) {
2079 if (CalcDerivatives(i,kTRUE)) return -1;
6526a72c 2080 m.derlocX[i] = fDerivativeLoc[i][0];
2081 m.derlocZ[i] = fDerivativeLoc[i][2];
2082 if (zeroX) zeroX = fDerivativeLoc[i][0]==0;
2083 if (zeroZ) zeroZ = fDerivativeLoc[i][2]==0;
7b85e477 2084 }
2085 // for (Int_t i=0; i<fNLocal; i++) AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
7b85e477 2086 //
6526a72c 2087 if (zeroX) {AliInfo("Aborting: zero local X derivatives!"); return -1;}
2088 if (zeroZ) {AliInfo("Aborting: zero local Z derivatives!"); return -1;}
2089 //
2090 int ifill = 0;
7b85e477 2091 //
2092 AliITSAlignMille2Module* endModule = fCurrentModule;
2093 //
6526a72c 2094 zeroX = zeroZ = kTRUE;
2095 Bool_t dfDone[kNParCh];
2096 for (int i=kNParCh;i--;) dfDone[i] = kFALSE;
2097 m.nModFilled = 0;
2098 //
2099 // special block for SDD derivatives
2100 Double_t jacobian[kNParChGeom];
2101 Int_t nmodTested = 0;
2102 //
7b85e477 2103 do {
6526a72c 2104 if (fCurrentModule->GetNParFree()==0) continue;
2105 nmodTested++;
2106 for (Int_t i=0; i<kNParChGeom; i++) { // common for all sensors: derivatives over geom params
2107 //
2108 if (!fUseGlobalDelta) dfDone[i] = kFALSE; // for global deltas the derivatives at diff. levels are different
2109 if (fCurrentModule->GetParOffset(i)<0) continue; // this parameter is not explicitly fitted
2110 if (!dfDone[i]) {
2111 if (CalcDerivatives(i,kFALSE)) return -1;
2112 else {
2113 dfDone[i] = kTRUE;
2114 if (zeroX) zeroX = fDerivativeGlo[i][0]==0;
2115 if (zeroZ) zeroZ = fDerivativeGlo[i][2]==0;
2116 }
2117 }
2118 //
2119 m.dergloX[ifill] = fDerivativeGlo[i][0];
2120 m.dergloZ[ifill] = fDerivativeGlo[i][2];
2121 m.parMilleID[ifill++] = fCurrentModule->GetParOffset(i);
2122 }
2123 //
2124 // specific for special sensors
2125 if ( fCurrentModule->IsSDD() &&
2126 (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0 ||
2127 fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) ) {
2128 //
2129 // assume for sensor local xloc = xloc0 + V0*dT0+dV*(T-T0)
2130 // where V0 and T are the nominal drift velocity, time and time0
2131 // and the dT0 and dV are the corrections:
2132 // dX/dT0 = dX/dxloc * dxloc/dT0 = dX/dxloc * V0
2133 // dX/dV = dX/dxloc * dxloc/dV = dX/dxloc * (T-T0)
2134 // IMPORTANT: the geom derivatives are over the SENSOR LOCAL parameters
2135 //
2136 if (!dfDone[AliITSAlignMille2Module::kDOFT0] || !dfDone[AliITSAlignMille2Module::kDOFDV]) {
2137 //
2138 double dXdxlocsens=0., dZdxlocsens=0.;
2139 //
2140 // if the current module is the sensor itself and we work with local params, then
2141 // we can directly take dX/dxloc_sens dZ/dxloc_sens
2142 if (!fUseGlobalDelta && fCurrentModule->GetVolumeID()==fCluster.GetVolumeID()) {
2143 if (dfDone[AliITSAlignMille2Module::kDOFTX]) {
2144 CalcDerivatives(AliITSAlignMille2Module::kDOFTX,kFALSE);
2145 dfDone[AliITSAlignMille2Module::kDOFTX] = kTRUE;
2146 }
2147 dXdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][0];
2148 dZdxlocsens = fDerivativeGlo[AliITSAlignMille2Module::kDOFTX][2];
2149 }
2150 //
2151 else { // need to perform some transformations
2152 // fetch the jacobian of the transformation from the sensors local frame to the frame
2153 // where the parameters are defined:
2154 // Global: dX/dxloc_sens = dX/dxgl*dxgl/dxloc_sens + ...dX/dphigl*dphigl/dxloc_sens
2155 if (fUseGlobalDelta) fCurrentModule->CalcDerivGloLoc(fCluster.GetVolumeID(),
2156 AliITSAlignMille2Module::kDOFTX, jacobian);
2157 // Local: dX/dxloc_sens = dX/dxcurr*dxcurr/dxloc_sens +..+dX/dphicurr * dphicurr/dxloc_sens
2158 else fCurrentModule->CalcDerivCurLoc(fCluster.GetVolumeID(),
2159 AliITSAlignMille2Module::kDOFTX, jacobian);
2160 //
2161 for (int j=0;j<kNParChGeom;j++) {
2162 // need global derivative even if the j-th param is locked
2163 if (!dfDone[j]) {CalcDerivatives(j,kFALSE); dfDone[j] = kTRUE;}
2164 dXdxlocsens += fDerivativeGlo[j][0] * jacobian[j];
2165 dZdxlocsens += fDerivativeGlo[j][2] * jacobian[j];
2166 }
2167 }
2168 //
2169 if (zeroX) zeroX = dXdxlocsens == 0;
2170 if (zeroZ) zeroZ = dZdxlocsens == 0;
2171 //
2172 double vdrift = GetVDriftSDD();
2173 double tdrift = GetTDriftSDD();
2174 //
2175 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0] = dXdxlocsens*vdrift;
2176 fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2] = dZdxlocsens*vdrift;
2177 dfDone[AliITSAlignMille2Module::kDOFT0] = kTRUE;
2178 //
2179 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0] = -dXdxlocsens*TMath::Sign(tdrift,vdrift);
2180 fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2] = -dZdxlocsens*TMath::Sign(tdrift,vdrift);
2181 dfDone[AliITSAlignMille2Module::kDOFDV] = kTRUE;
2182 //
2183 }
2184 //
2185 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0)>=0) {
2186 m.dergloX[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][0];
2187 m.dergloZ[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFT0][2];
2188 m.parMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFT0);
2189 }
2190 //
2191 if (fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV)>=0) {
2192 m.dergloX[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][0];
2193 m.dergloZ[ifill] = fDerivativeGlo[AliITSAlignMille2Module::kDOFDV][2];
2194 m.parMilleID[ifill++] = fCurrentModule->GetParOffset(AliITSAlignMille2Module::kDOFDV);
7b85e477 2195 }
2196 }
7b85e477 2197 //
6526a72c 2198 m.moduleID[m.nModFilled++] = fCurrentModule->GetUniqueID();
7b85e477 2199 } while( (fCurrentModule=fCurrentModule->GetParent()) );
2200 //
6526a72c 2201 if (nmodTested>0 && zeroX) {AliInfo("Aborting: zero global X derivatives!");return -1;}
2202 if (nmodTested>0 && zeroZ) {AliInfo("Aborting: zero global Z derivatives!");return -1;}
2203 //
7b85e477 2204 // ok, can copy to m
2205 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
2206 m.measX = fMeasLoc[0]-fPintLoc0[0];
2207 m.sigmaX = fSigmaLoc[0];
2208 //
2209 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
2210 m.measZ = fMeasLoc[2]-fPintLoc0[2];
2211 m.sigmaZ = fSigmaLoc[2];
2212 //
6526a72c 2213 m.nGlobFilled = ifill;
7b85e477 2214 fCurrentModule = endModule;
2215 //
6526a72c 2216 return Int_t(!zeroX && !zeroZ);
7b85e477 2217}
2218
6526a72c 2219//________________________________________________________________________________________________________
2220void AliITSAlignMille2::SetLocalEquations(const Mille2Data *marr, Int_t neq)
2221{
7b85e477 2222 /// Set local equations with data stored in m
2223 /// return 0 if success
2224 //
2225 for (Int_t j=0; j<neq; j++) {
2226 //
2227 const Mille2Data &m = marr[j];
2228 //
2229 // set equation for Xloc coordinate
2230 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",m.measX, m.sigmaX));
2231 for (int i=fNLocal; i--;) SetLocalDerivative( i, m.derlocX[i] );
6526a72c 2232 for (int i=m.nGlobFilled;i--;) SetGlobalDerivative( m.parMilleID[i] , m.dergloX[i] );
7b85e477 2233 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.measX, m.sigmaX);
2234 //
2235 // set equation for Zloc coordinate
2236 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",m.measZ, m.sigmaZ));
2237 for (int i=fNLocal; i--;) SetLocalDerivative( i, m.derlocZ[i] );
6526a72c 2238 for (int i=m.nGlobFilled;i--;) SetGlobalDerivative( m.parMilleID[i] , m.dergloZ[i] );
7b85e477 2239 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m.measZ, m.sigmaZ);
6526a72c 2240 //
2241 for (int i=m.nModFilled;i--;) GetMilleModule(m.moduleID[i])->IncNProcessedPoints();
2242 //
7b85e477 2243 }
2244}
2245
6526a72c 2246//________________________________________________________________________________________________________
2247Int_t AliITSAlignMille2::GlobalFit()
2248{
7b85e477 2249 /// Call global fit; Global parameters are stored in parameters
6526a72c 2250 if (!fIsMilleInit) Init();
2251 //
2252 ApplyPreConstraints();
2253 int res = fMillepede->GlobalFit();
7b85e477 2254 AliInfo(Form("%s fitting global parameters!",res ? "Done":"Failed"));
6526a72c 2255 if (res) {
2256 // fetch the parameters
2257 for (int imd=fNModules;imd--;) {
2258 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2259 int nprocp = 0;
2260 for (int ip=mod->GetNParTot();ip--;) {
2261 int idp = mod->GetParOffset(ip);
2262 if (idp<0) continue; // was not in the explicit fit
2263 mod->SetParVal(ip,fMillepede->GetFinalParam(idp));
2264 mod->SetParErr(ip,fMillepede->GetFinalError(idp));
2265 int np = fMillepede->GetProcessedPoints(idp);
2266 if (TMath::Abs(np)>TMath::Abs(nprocp)) nprocp = np;
2267 }
2268 if (!mod->GetNProcessedPoints()) mod->SetNProcessedPoints(nprocp);
2269 }
7b85e477 2270
7b85e477 2271 }
6526a72c 2272 ApplyPostConstraints();
2273 return res;
7b85e477 2274}
2275
6526a72c 2276//________________________________________________________________________________________________________
2277void AliITSAlignMille2::PrintGlobalParameters()
2278{
7b85e477 2279 /// Print global parameters
2280 if (!fIsMilleInit) {
2281 AliInfo("Millepede has not been initialized!");
2282 return;
2283 }
2284 fMillepede->PrintGlobalParameters();
2285}
2286
6526a72c 2287//________________________________________________________________________________________________________
7b85e477 2288Int_t AliITSAlignMille2::LoadSuperModuleFile(const Char_t *sfile)
2289{
2290 // load definitions of supermodules from a root file
2291 // return 0 if success
2292
2293 TFile *smf=TFile::Open(sfile);
2294 if (!smf->IsOpen()) {
2295 AliInfo(Form("Cannot open supermodule file %s",sfile));
2296 return -1;
2297 }
2298
2299 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
2300 if (!sma) {
2301 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
2302 return -2;
2303 }
2304 Int_t nsma=sma->GetEntriesFast();
2305 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
2306 //
2307 Char_t st[250];
2308 char symname[150];
2309 UShort_t volid;
2310 TGeoHMatrix m;
2311 //
2312 for (Int_t i=0; i<nsma; i++) {
2313 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
2314 volid=a->GetVolUID();
2315 strcpy(st,a->GetSymName());
2316 a->GetMatrix(m);
2317 //
2318 sscanf(st,"%s",symname);
6526a72c 2319 //
7b85e477 2320 // decode module list
2321 char *stp=strstr(st,"ModuleList:");
2322 if (!stp) return -3;
2323 stp += 11;
2324 int idx[2200];
2325 char spp[200]; int jp=0;
2326 char cl[20];
2327 strcpy(st,stp);
2328 int l=strlen(st);
2329 int j=0;
2330 int n=0;
2331 //
2332 while (j<=l) {
2333 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
2334 spp[jp]=0;
2335 jp=0;
2336 if (strlen(spp)) {
2337 int k=strcspn(spp,"-");
2338 if (k<int(strlen(spp))) { // c'e' il -
2339 strcpy(cl,&(spp[k+1]));
2340 spp[k]=0;
2341 int ifrom=atoi(spp); int ito=atoi(cl);
2342 for (int b=ifrom; b<=ito; b++) {
2343 idx[n]=b;
2344 n++;
2345 }
2346 }
2347 else { // numerillo singolo
2348 idx[n]=atoi(spp);
2349 n++;
2350 }
2351 }
2352 }
2353 else {
2354 spp[jp]=st[j];
2355 jp++;
2356 }
2357 j++;
2358 }
2359 UShort_t volidsv[2198];
2360 for (j=0;j<n;j++) {
2361 volidsv[j]=AliITSAlignMille2Module::GetVolumeIDFromIndex(idx[j]);
2362 if (!volidsv[j]) {
2363 AliInfo(Form("Index %d not valid (range 0->%d)",idx[j],kMaxITSSensID));
2364 return -5;
2365 }
2366 }
2367 Int_t smindex=int(2198+volid-14336); // virtual index
6526a72c 2368 //
7b85e477 2369 fSuperModule.AddAtAndExpand(new AliITSAlignMille2Module(smindex,volid,symname,&m,n,volidsv),fNSuperModules);
2370 //
2371 fNSuperModules++;
2372 }
6526a72c 2373 //
7b85e477 2374 smf->Close();
2375 //
2376 return 0;
2377}
2378
6526a72c 2379//________________________________________________________________________________________________________
2380void AliITSAlignMille2::ConstrainModuleSubUnitsMean(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2381{
2382 // require that sum of modifications for the childs of this module is = val, i.e.
2383 // the internal corrections moves the module as a whole by fixed value (0 by default).
2384 // pattern is the bit pattern for the parameters to constrain
2385 //
6526a72c 2386 if (fIsMilleInit) {
2387 AliInfo("Millepede has been already initialized: no constrain may be added!");
2388 return;
7b85e477 2389 }
6526a72c 2390 if (!GetMilleModule(idm)->GetNChildren()) return;
2391 TString nm = "cstrSUMean";
2392 nm += GetNConstraints();
2393 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2394 idm,val,pattern);
2395 cstr->SetConstraintID(GetNConstraints());
2396 fConstraints.Add(cstr);
7b85e477 2397}
2398
6526a72c 2399//________________________________________________________________________________________________________
2400void AliITSAlignMille2::ConstrainModuleSubUnitsMedian(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2401{
2402 // require that median of the modifications for the childs of this module is = val, i.e.
2403 // the internal corrections moves the module as a whole by fixed value (0 by default)
2404 // module the outliers.
2405 // pattern is the bit pattern for the parameters to constrain
2406 // The difference between the mean and the median will be transfered to the parent
6526a72c 2407 if (fIsMilleInit) {
2408 AliInfo("Millepede has been already initialized: no constrain may be added!");
2409 return;
2410 }
2411 if (!GetMilleModule(idm)->GetNChildren()) return;
2412 TString nm = "cstrSUMed";
2413 nm += GetNConstraints();
2414 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2415 idm,val,pattern);
2416 cstr->SetConstraintID(GetNConstraints());
2417 fConstraints.Add(cstr);
2418}
2419
2420//________________________________________________________________________________________________________
2421void AliITSAlignMille2::ConstrainOrphansMean(Double_t val, UInt_t pattern)
2422{
2423 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2424 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2425 // pattern is the bit pattern for the parameters to constrain
7b85e477 2426 //
6526a72c 2427 if (fIsMilleInit) {
2428 AliInfo("Millepede has been already initialized: no constrain may be added!");
2429 return;
2430 }
2431 TString nm = "cstrOMean";
2432 nm += GetNConstraints();
2433 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMean,
2434 -1,val,pattern);
2435 cstr->SetConstraintID(GetNConstraints());
2436 fConstraints.Add(cstr);
2437}
2438
2439//________________________________________________________________________________________________________
2440void AliITSAlignMille2::ConstrainOrphansMedian(Double_t val, UInt_t pattern)
2441{
2442 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2443 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2444 // pattern is the bit pattern for the parameters to constrain
7b85e477 2445 //
6526a72c 2446 if (fIsMilleInit) {
2447 AliInfo("Millepede has been already initialized: no constrain may be added!");
2448 return;
7b85e477 2449 }
6526a72c 2450 TString nm = "cstrOMed";
2451 nm += GetNConstraints();
2452 AliITSAlignMille2Constraint *cstr = new AliITSAlignMille2Constraint(nm.Data(),AliITSAlignMille2Constraint::kTypeMedian,
2453 -1,val,pattern);
2454 cstr->SetConstraintID(GetNConstraints());
2455 fConstraints.Add(cstr);
2456}
2457
2458//________________________________________________________________________________________________________
2459void AliITSAlignMille2::ConstrainLocal(const Char_t* name,Double_t *parcf,Int_t npar,Double_t val,Double_t err)
2460{
2461 if (fIsMilleInit) {
2462 AliInfo("Millepede has been already initialized: no constrain may be added!");
2463 return;
2464 }
2465 AliITSAlignMille2ConstrArray *cstr = new AliITSAlignMille2ConstrArray(name,parcf,npar,val,err);
2466 cstr->SetConstraintID(GetNConstraints());
2467 fConstraints.Add(cstr);
2468}
2469
2470//________________________________________________________________________________________________________
2471void AliITSAlignMille2::ApplyGaussianConstraint(AliITSAlignMille2ConstrArray* cstr)
2472{
2473 // apply the constraint on the local corrections of a list of modules
2474 int nmod = cstr->GetNModules();
2475 double jacobian[AliITSAlignMille2Module::kMaxParGeom][AliITSAlignMille2Module::kMaxParGeom];
7b85e477 2476 //
6526a72c 2477 for (int imd=nmod;imd--;) {
2478 int modID = cstr->GetModuleID(imd);
2479 AliITSAlignMille2Module* mod = GetMilleModule(modID);
2480 ResetLocalEquation();
2481 int nadded = 0;
2482 double value = cstr->GetValue();
2483 double sigma = cstr->GetError();
2484 //
2485 // in case the reference (survey) deltas were imposed for Gaussian constraints
2486 // already accumulated corrections: they must be subtracted from the constraint value.
2487 if (IsConstraintWrtRef()) {
2488 //
2489 Double_t precal[AliITSAlignMille2Module::kMaxParTot];
2490 Double_t refcal[AliITSAlignMille2Module::kMaxParTot];
2491 for (int ip=AliITSAlignMille2Module::kMaxParTot;ip--;) {precal[ip]=0; refcal[ip] = 0.;}
2492 //
2493 // check if there was a reference delta provided for this module
2494 AliAlignObjParams* parref = GetConstrRefObject(mod->GetName());
2495 if (parref) parref->GetPars(refcal, refcal+3); // found reference delta
2496 //
2497 // extract already applied local corrections for this module
2498 if (fPrealignment) {
2499 //
2500 AliAlignObjParams *preo = GetPrealignedObject(mod->GetName());
2501 if (preo) {
2502 TGeoHMatrix preMat,tmpMat = *mod->GetMatrix(); // Delta_Glob * Delta_Glob_Par * M
2503 preo->GetMatrix(preMat); // Delta_Glob
2504 preMat.MultiplyLeft( &tmpMat.Inverse() ); // M^-1 * Delta_Glob_Par^-1 = (Delta_Glob_Par * M)^-1
2505 tmpMat.MultiplyLeft( &preMat ); // (Delta_Glob_Par * M)^-1 * Delta_Glob * Delta_Glob_Par * M = Delta_loc
2506 AliAlignObjParams algob;
2507 algob.SetMatrix(tmpMat);
2508 algob.GetPars(precal,precal+3); // local corrections for geometry
2509 }
2510 }
2511 //
2512 // subtract the contribution to constraint from precalibration
2513 for (int ipar=cstr->GetNCoeffs();ipar--;) value += (refcal[ipar]-precal[ipar])*cstr->GetCoeff(ipar);
2514 //
2515 }
2516 //
2517 if (fUseGlobalDelta) mod->CalcDerivLocGlo(&jacobian[0][0]);
2518 //
2519 for (int ipar=cstr->GetNCoeffs();ipar--;) {
2520 double coef = cstr->GetCoeff(ipar);
2521 if (coef==0) continue;
2522 //
2523 if (!fUseGlobalDelta || ipar>= AliITSAlignMille2Module::kMaxParGeom) { //
2524 // we are working with local params or if the given param is not related to geometry,
2525 // apply the constraint directly
2526 int parPos = mod->GetParOffset(ipar);
2527 if (parPos<0) continue; // not in the fit
2528 fGlobalDerivatives[parPos] += coef;
2529 nadded++;
2530 }
2531 else { // we are working with global params, while the constraint is on local ones -> jacobian
2532 for (int jpar=AliITSAlignMille2Module::kMaxParGeom;jpar--;) {
2533 int parPos = mod->GetParOffset(jpar);
2534 if (parPos<0) continue;
2535 fGlobalDerivatives[parPos] += coef*jacobian[ipar][jpar];
2536 nadded++;
2537 }
2538 }
2539 }
2540 if (nadded) AddConstraint(fGlobalDerivatives, value, sigma);
2541 }
7b85e477 2542 //
6526a72c 2543}
2544
2545//________________________________________________________________________________________________________
2546void AliITSAlignMille2::ApplyPreConstraints()
2547{
2548 int nconstr = GetNConstraints();
2549 for (int i=0;i<nconstr;i++) {
2550 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
7b85e477 2551 //
6526a72c 2552 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) {
2553 ApplyGaussianConstraint( (AliITSAlignMille2ConstrArray*)cstr);
2554 continue;
2555 }
7b85e477 2556 //
6526a72c 2557 if (cstr->GetType() == AliITSAlignMille2Constraint::kTypeMedian) continue; // post type constraint
2558 //
2559 if (!fUseGlobalDelta) continue; // mean/med constraints must be applied to global deltas
2560 // apply constraint on the mean's before the fit
2561 int imd = cstr->GetModuleID();
2562 if (imd>=0) {
2563 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2564 UInt_t pattern = 0;
2565 for (int ipar=mod->GetNParTot();ipar--;) {
2566 if (!cstr->IncludesParam(ipar)) continue;
2567 if (mod->GetParOffset(ipar)<0) continue; // parameter is not in the explicit fit -> post constraint
2568 pattern |= 0x1<<ipar;
2569 cstr->SetApplied(ipar);
2570 }
2571 ConstrainModuleSubUnits(imd,cstr->GetValue(),pattern);
2572 //
2573 }
2574 else if (!PseudoParentsAllowed()) {
2575 ConstrainOrphans(cstr->GetValue(),(UInt_t)cstr->GetPattern());
2576 cstr->SetApplied(-1);
2577 }
7b85e477 2578 }
6526a72c 2579}
2580
2581//________________________________________________________________________________________________________
2582void AliITSAlignMille2::ApplyPostConstraints()
2583{
2584 int nconstr = GetNConstraints();
2585 Bool_t convGlo = kFALSE;
2586 // check if there is something to do
2587 int ntodo = 0;
2588 for (int i=0;i<nconstr;i++) {
2589 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2590 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2591 if (cstr->GetRemainingPattern() == 0) continue;
2592 ntodo++;
2593 }
2594 if (!ntodo) return;
7b85e477 2595 //
6526a72c 2596 if (!fUseGlobalDelta) { // need to convert to global params
2597 ConvertParamsToGlobal();
2598 convGlo = kTRUE;
2599 }
7b85e477 2600 //
6526a72c 2601 for (int i=0;i<nconstr;i++) {
2602 AliITSAlignMille2Constraint* cstr = GetConstraint(i);
2603 if (cstr->GetType() == AliITSAlignMille2ConstrArray::kTypeGaussian) continue;
2604 //
2605 int imd = cstr->GetModuleID();
2606 //
2607 if (imd>=0) {
2608 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2609 UInt_t pattern = 0;
2610 for (int ipar=mod->GetNParTot();ipar--;) {
2611 if (cstr->IsApplied(ipar)) continue;
2612 if (!cstr->IncludesParam(ipar)) continue;
2613 if (!mod->IsFreeDOF(ipar)) continue; // parameter is fixed, will not apply constraint
2614 pattern |= 0x1<<ipar;
2615 cstr->SetApplied(ipar);
2616 }
2617 if (pattern) PostConstrainModuleSubUnits(cstr->GetType(),cstr->GetModuleID(),cstr->GetValue(),pattern);
2618 //
2619 }
2620 else if (PseudoParentsAllowed()) {
2621 UInt_t pattern = (UInt_t)cstr->GetRemainingPattern();
2622 PostConstrainOrphans(cstr->GetType(),cstr->GetValue(),pattern);
2623 cstr->SetApplied(-1);
2624 }
2625 }
2626 // if there was a conversion, rewind it
2627 if (convGlo) ConvertParamsToLocal();
2628 //
7b85e477 2629}
2630
6526a72c 2631//________________________________________________________________________________________________________
2632void AliITSAlignMille2::ConstrainModuleSubUnits(Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2633{
6526a72c 2634 // require that sum of modifications for the childs of this module is = val, i.e.
2635 // the internal corrections moves the module as a whole by fixed value (0 by default).
7b85e477 2636 // pattern is the bit pattern for the parameters to constrain
2637 //
7b85e477 2638 //
6526a72c 2639 AliITSAlignMille2Module* mod = GetMilleModule(idm);
7b85e477 2640 //
6526a72c 2641 for (int ip=0;ip<kNParCh;ip++) {
2642 if ( !((pattern>>ip)&0x1) /*|| !parent->IsFreeDOF(ip)*/) continue;
2643 ResetLocalEquation();
2644 int nadd = 0;
2645 for (int ich=mod->GetNChildren();ich--;) {
2646 int idpar = ((AliITSAlignMille2Module*)mod->GetChild(ich))->GetParOffset(ip);
2647 if (idpar<0) continue;
2648 fGlobalDerivatives[idpar] = 1.0;
2649 nadd++;
2650 }
2651 //
2652 if (nadd>0) {
2653 AddConstraint(fGlobalDerivatives,val);
2654 AliInfo(Form("Constrained param %d for %d submodules of module #%d: %s",ip,nadd,idm,mod->GetName()));
2655 }
7b85e477 2656 }
7b85e477 2657 //
6526a72c 2658}
2659
2660//________________________________________________________________________________________________________
2661void AliITSAlignMille2::ConstrainOrphans(Double_t val, UInt_t pattern)
2662{
2663 // require that median of the modifications for the supermodules which have no parents is = val, i.e.
2664 // the corrections moves the whole setup by fixed value (0 by default) modulo the outliers.
2665 // pattern is the bit pattern for the parameters to constrain
2666 //
7b85e477 2667 for (int ip=0;ip<kNParCh;ip++) {
6526a72c 2668 //
2669 if ( !((pattern>>ip)&0x1) ) continue;
7b85e477 2670 ResetLocalEquation();
6526a72c 2671 int nadd = 0;
2672 for (int imd=fNModules;imd--;) {
2673 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2674 if (mod->GetParent()) continue; // this is not an orphan
2675 int idpar = mod->GetParOffset(ip);
2676 if (idpar<0) continue;
2677 fGlobalDerivatives[idpar] = 1.0;
2678 nadd++;
2679 }
2680 if (nadd>0) {
2681 AddConstraint(fGlobalDerivatives,val);
2682 AliInfo(Form("Constrained param %d for %d orphan modules",ip,nadd));
7b85e477 2683 }
7b85e477 2684 }
2685 //
7b85e477 2686 //
2687}
6526a72c 2688
2689//________________________________________________________________________________________________________
2690void AliITSAlignMille2::PostConstrainModuleSubUnits(Int_t type,Int_t idm, Double_t val, UInt_t pattern)
7b85e477 2691{
6526a72c 2692 // require that median or mean of the modifications for the childs of this module is = val, i.e.
2693 // the internal corrections moves the module as a whole by fixed value (0 by default)
2694 // module the outliers.
7b85e477 2695 // pattern is the bit pattern for the parameters to constrain
6526a72c 2696 // The difference between the mean and the median will be transfered to the parent
7b85e477 2697 //
6526a72c 2698 AliITSAlignMille2Module* parent = GetMilleModule(idm);
2699 int nc = parent->GetNChildren();
7b85e477 2700 //
6526a72c 2701 double *tmpArr = new double[nc];
7b85e477 2702 //
2703 for (int ip=0;ip<kNParCh;ip++) {
6526a72c 2704 int npc = 0;
2705 if ( !((pattern>>ip)&0x1) || !parent->IsFreeDOF(ip)) continue;
2706 // compute the mean and median of the deltas
2707 int nfree = 0;
2708 for (int ich=nc;ich--;) {
2709 AliITSAlignMille2Module* child = parent->GetChild(ich);
2710 // if (!child->IsFreeDOF(ip)) continue;
2711 tmpArr[nfree++] = child->GetParVal(ip);
2712 }
2713 double median=0,mean=0;
2714 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
2715 mean += tmpArr[ic0];
2716 for (int ic1=ic0+1;ic1<nfree;ic1++)
2717 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
2718 }
7b85e477 2719 //
6526a72c 2720 int kmed = nfree/2;
2721 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
2722 if (nfree>0) mean /= nfree;
2723 //
2724 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
2725 //
2726 for (int ich=nc;ich--;) {
2727 AliITSAlignMille2Module* child = parent->GetChild(ich);
2728 // if (!child->IsFreeDOF(ip)) continue;
2729 child->SetParVal(ip, child->GetParVal(ip) + shift);
2730 npc++;
7b85e477 2731 }
6526a72c 2732 //
2733 parent->SetParVal(ip, parent->GetParVal(ip) - shift);
2734 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d children of module %d: %s",
2735 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
2736 ip,npc,idm,parent->GetName()));
7b85e477 2737 }
6526a72c 2738 delete[] tmpArr;
7b85e477 2739 //
7b85e477 2740 //
2741}
2742
6526a72c 2743//________________________________________________________________________________________________________
2744void AliITSAlignMille2::PostConstrainOrphans(Int_t type,Double_t val, UInt_t pattern)
7b85e477 2745{
6526a72c 2746 // require that median or mean of modifications for the supermodules which have no parents is = val, i.e.
2747 // the corrections moves the whole setup by fixed value (0 by default).
2748 // pattern is the bit pattern for the parameters to constrain
7b85e477 2749 //
6526a72c 2750 int nc = fNModules;
7b85e477 2751 //
6526a72c 2752 int norph = 0;
2753 for (int ich=nc;ich--;) if (!GetMilleModule(ich)->GetParent()) norph ++;
2754 if (!norph) return;
2755 double *tmpArr = new double[norph];
7b85e477 2756 //
6526a72c 2757 for (int ip=0;ip<kNParCh;ip++) {
2758 int npc = 0;
2759 if ( !((pattern>>ip)&0x1)) continue;
2760 // compute the mean and median of the deltas
2761 int nfree = 0;
2762 for (int ich=nc;ich--;) {
2763 AliITSAlignMille2Module* child = GetMilleModule(ich);
2764 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
2765 if (child->GetParent()) continue;
2766 tmpArr[nfree++] = child->GetParVal(ip);
2767 }
2768 double median=0,mean=0;
2769 for (int ic0=0;ic0<nfree;ic0++) {// order the deltas
2770 mean += tmpArr[ic0];
2771 for (int ic1=ic0+1;ic1<nfree;ic1++)
2772 if (tmpArr[ic0]>tmpArr[ic1]) {double tv=tmpArr[ic0]; tmpArr[ic0]=tmpArr[ic1]; tmpArr[ic1]=tv;}
2773 }
2774 //
2775 int kmed = nfree/2;
2776 median = (tmpArr[kmed]+tmpArr[nfree-kmed-1])/2.;
2777 if (nfree>0) mean /= nfree;
2778 //
2779 double shift = val - (type==AliITSAlignMille2Constraint::kTypeMean ? mean : median);
2780 //
2781 for (int ich=nc;ich--;) {
2782 AliITSAlignMille2Module* child = GetMilleModule(ich);
2783 // if (child->GetParent() || !child->IsFreeDOF(ip)) continue;
2784 if (child->GetParent()) continue;
2785 child->SetParVal(ip, child->GetParVal(ip) + shift);
2786 npc++;
7b85e477 2787 }
6526a72c 2788 //
2789 AliInfo(Form("%s constraint: added %f shift to param[%d] of %d orphan modules",
2790 type==AliITSAlignMille2Constraint::kTypeMean ? "MEAN" : "MEDIAN",shift,
2791 ip,npc));
7b85e477 2792 }
6526a72c 2793 delete[] tmpArr;
7b85e477 2794 //
6526a72c 2795}
2796
2797//________________________________________________________________________________________________________
2798Bool_t AliITSAlignMille2::IsParModConstrained(AliITSAlignMille2Module* mod,Int_t par, Bool_t &meanmed, Bool_t &gaussian) const
2799{
2800 // check if par of the module participates in some constraint, and set the flag for their types
2801 meanmed = gaussian = kFALSE;
2802 //
2803 if ( mod->IsParConstrained(par) ) gaussian = kTRUE; // direct constraint on this param
2804 //
2805 for (int icstr=GetNConstraints();icstr--;) {
2806 AliITSAlignMille2Constraint* cstr = GetConstraint(icstr);
2807 //
2808 if (!cstr->IncludesModPar(mod,par)) continue;
2809 if (cstr->GetType()==AliITSAlignMille2ConstrArray::kTypeGaussian) gaussian = kTRUE;
2810 else meanmed = kTRUE;
2811 //
2812 if (meanmed && gaussian) break; // no sense to check further
7b85e477 2813 }
6526a72c 2814 //
2815 return meanmed||gaussian;
2816}
2817
2818//________________________________________________________________________________________________________
2819Bool_t AliITSAlignMille2::IsParModFamilyVaried(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
2820{
2821 // check if parameter par is varied for this module or its children up to the level depth
2822 if (depth<0) return kFALSE;
2823 if (mod->GetParOffset(par)>=0) return kTRUE;
2824 for (int icld=mod->GetNChildren();icld--;) {
2825 AliITSAlignMille2Module* child = mod->GetChild(icld);
2826 if (IsParModFamilyVaried(child, par, depth-1)) return kTRUE;
7b85e477 2827 }
6526a72c 2828 return kFALSE;
7b85e477 2829 //
6526a72c 2830}
2831
2832/*
2833//________________________________________________________________________________________________________
2834Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
2835{
2836 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
2837 if (depth<0) return kTRUE;
2838 for (int icld=mod->GetNChildren();icld--;) {
2839 AliITSAlignMille2Module* child = mod->GetChild(icld);
2840 //if (child->GetParOffset(par)<0) continue; // fixed
2841 Bool_t cstMM=kFALSE,cstGS=kFALSE;
2842 // does this child have gaussian constraint ?
2843 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
2844 // check its children
2845 if (!IsParFamilyFree(child,par,depth-1)) return kTRUE;
2846 }
2847 return kFALSE;
2848 //
2849}
2850*/
2851
2852//________________________________________________________________________________________________________
2853Bool_t AliITSAlignMille2::IsParFamilyFree(AliITSAlignMille2Module* mod,Int_t par,Int_t depth) const
2854{
2855 // check if parameter par is varied and is not subjected to gaussian constraint for the children up to the level depth
2856 if (depth<0) return kFALSE;
2857 for (int icld=mod->GetNChildren();icld--;) {
2858 AliITSAlignMille2Module* child = mod->GetChild(icld);
2859 //if (child->GetParOffset(par)<0) continue; // fixed
2860 Bool_t cstMM=kFALSE,cstGS=kFALSE;
2861 // does this child have gaussian constraint ?
2862 if (!IsParModConstrained(child,par,cstMM,cstGS) || !cstGS ) return kTRUE;
2863 // check its children
2864 if (IsParFamilyFree(child,par,depth-1)) return kTRUE;
2865 }
2866 return kFALSE;
7b85e477 2867 //
2868}
6526a72c 2869
2870//________________________________________________________________________________________________________
2871Double_t AliITSAlignMille2::GetTDriftSDD() const
2872{
2873 double t = fCluster.GetDriftTime();
2874 return t - fDriftTime0[ fCluster.GetUniqueID() ];
2875}
2876
2877//________________________________________________________________________________________________________
2878Double_t AliITSAlignMille2::GetVDriftSDD() const
2879{
2880 return fDriftSpeed[ fCluster.GetUniqueID() ];
2881}
2882
2883//________________________________________________________________________________________________________
2884Bool_t AliITSAlignMille2::FixedOrphans() const
2885{
2886 // are there fixed modules with no parent (normally in such a case
2887 // the constraints on the orphans should not be applied
2888 if (!IsConfigured()) {
2889 AliInfo("Still not configured");
2890 return kFALSE;
2891 }
2892 for (int i=0;i<fNModules;i++) {
2893 AliITSAlignMille2Module* md = GetMilleModule(i);
2894 if (md->GetParent()==0 && md->GetNParFree()==0) return kTRUE;
2895 }
2896 return kFALSE;
2897}
2898
2899//________________________________________________________________________________________________________
2900void AliITSAlignMille2::ConvertParamsToGlobal()
2901{
2902 double pars[AliITSAlignMille2Module::kMaxParGeom];
2903 for (int imd=fNModules;imd--;) {
2904 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2905 if (mod->GeomParamsGlobal()) continue;
2906 mod->GetGeomParamsGlo(pars);
2907 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
2908 mod->SetGeomParamsGlobal(kTRUE);
2909 }
2910}
2911
2912//________________________________________________________________________________________________________
2913void AliITSAlignMille2::ConvertParamsToLocal()
2914{
2915 double pars[AliITSAlignMille2Module::kMaxParGeom];
2916 for (int imd=fNModules;imd--;) {
2917 AliITSAlignMille2Module* mod = GetMilleModule(imd);
2918 if (!mod->GeomParamsGlobal()) continue;
2919 mod->GetGeomParamsLoc(pars);
2920 mod->SetParVals(pars,AliITSAlignMille2Module::kMaxParGeom);
2921 mod->SetGeomParamsGlobal(kFALSE);
2922 }
2923}
2924
2925