Updated macro
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille.cxx
CommitLineData
483fd7d8 1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//-----------------------------------------------------------------------------
19/// \class AliITSAlignMille
20/// Alignment class fro 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)
28//-----------------------------------------------------------------------------
29
30#include <TF1.h>
0856765c 31#include <TFile.h>
32#include <TClonesArray.h>
483fd7d8 33#include <TGraph.h>
34#include <TGeoMatrix.h>
35#include <TMath.h>
0856765c 36#include <TGraphErrors.h>
483fd7d8 37
0856765c 38#include "AliITSAlignMilleModule.h"
483fd7d8 39#include "AliITSAlignMille.h"
40#include "AliITSgeomTGeo.h"
41#include "AliGeomManager.h"
42#include "AliMillepede.h"
43#include "AliTrackPointArray.h"
44#include "AliAlignObjParams.h"
45#include "AliLog.h"
46#include "TSystem.h" // come si fa?
47
48/// \cond CLASSIMP
49ClassImp(AliITSAlignMille)
50/// \endcond
51
52Int_t AliITSAlignMille::fgNDetElem = ITSMILLE_NDETELEM;
53Int_t AliITSAlignMille::fgNParCh = ITSMILLE_NPARCH;
54
55AliITSAlignMille::AliITSAlignMille(const Char_t *configFilename, Bool_t initmille)
56 : TObject(),
57 fMillepede(0),
58 fStartFac(16.),
59 fResCutInitial(100.),
60 fResCut(100.),
61 fNGlobal(ITSMILLE_NDETELEM*ITSMILLE_NPARCH),
62 fNLocal(ITSMILLE_NLOCAL),
63 fNStdDev(ITSMILLE_NSTDEV),
64 fIsMilleInit(kFALSE),
65 fParSigTranslations(0.0100),
66 fParSigRotations(0.1),
67 fTrack(NULL),
68 fCluster(),
69 fGlobalDerivatives(NULL),
483fd7d8 70 fTempAlignObj(NULL),
71 fDerivativeXLoc(0),
72 fDerivativeZLoc(0),
73 fDeltaPar(0),
74 fMinNPtsPerTrack(3),
0856765c 75 fInitTrackParamsMeth(1),
483fd7d8 76 fGeometryFileName("geometry.root"),
0856765c 77 fPreAlignmentFileName(""),
483fd7d8 78 fGeoManager(0),
79 fCurrentModuleIndex(0),
80 fCurrentModuleInternalIndex(0),
0856765c 81 fCurrentSensVolIndex(0),
483fd7d8 82 fNModules(0),
83 fUseLocalShifts(kTRUE),
0856765c 84 fUseSuperModules(kFALSE),
85 fUsePreAlignment(kFALSE),
86 fNSuperModules(0),
483fd7d8 87 fCurrentModuleHMatrix(NULL)
88{
89 /// main constructor that takes input from configuration file
90
91 fMillepede = new AliMillepede();
92 fGlobalDerivatives = new Double_t[fNGlobal];
0856765c 93 //fTempHMat = new TGeoHMatrix;
94 //fCurrentModuleHMatrix = new TGeoHMatrix;
483fd7d8 95
96 Int_t lc=LoadConfig(configFilename);
97 if (lc) {
98 AliInfo(Form("Error %d loading configuration from %s",lc,configFilename));
99 }
100 else {
101 if (initmille && fNGlobal<20000) {
102 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
103 Init(fNGlobal, fNLocal, fNStdDev);
104 ResetLocalEquation();
105 AliInfo("Parameters initialized to zero");
106 }
107 else {
108 AliInfo("Millepede has not been initialized ... ");
109 }
110 }
111
112 fDeltaPar=0.0; // not used at the moment - to be checked later...
113
114}
115
116AliITSAlignMille::~AliITSAlignMille() {
117 /// Destructor
118 if (fMillepede) delete fMillepede;
119 delete [] fGlobalDerivatives;
0856765c 120 //delete fCurrentModuleHMatrix;
121 //delete fTempHMat;
122 for (int i=0; i<fNModules; i++) delete fMilleModule[i];
123 for (int i=0; i<fNSuperModules; i++) delete fSuperModule[i];
483fd7d8 124}
125
126///////////////////////////////////////////////////////////////////////
127
128Int_t AliITSAlignMille::LoadConfig(const Char_t *cfile) {
129 /// return 0 if success
130 /// 1 if error in module index or voluid
131
132 FILE *pfc=fopen(cfile,"r");
133 if (!pfc) return -1;
134
135 Char_t st[200],st2[200];
136 Char_t tmp[100];
137 Int_t idx,itx,ity,itz,ith,ips,iph;
138 Float_t f1;
139 UShort_t voluid;
140 Int_t nmod=0;
141
142 while (fgets(st,200,pfc)) {
143
144 // skip comments
145 for (int i=0; i<int(strlen(st)); i++) {
146 if (st[i]=='#') st[i]=0;
147 }
148
149 if (strstr(st,"GEOMETRY_FILE")) {
150 sscanf(st,"%s %s",tmp,st2);
151 if (gSystem->AccessPathName(st2)) {
152 AliInfo("*** WARNING! *** geometry file not found! ");
153 return -1;
154 }
155 fGeometryFileName=st2;
156 InitGeometry();
157 }
158
0856765c 159 if (strstr(st,"PREALIGNMENT_FILE")) {
160 sscanf(st,"%s %s",tmp,st2);
161 if (gSystem->AccessPathName(st2)) {
162 AliInfo("*** WARNING! *** prealignment file not found! ");
163 return -1;
164 }
165 fPreAlignmentFileName=st2;
166 itx=ApplyToGeometry();
167 if (itx) {
168 AliInfo(Form("*** WARNING! *** error %d reading prealignment file! ",itx));
169 return -6;
170 }
171 }
172
173 if (strstr(st,"SUPERMODULE_FILE")) {
174 sscanf(st,"%s %s",tmp,st2);
175 if (gSystem->AccessPathName(st2)) {
176 AliInfo("*** WARNING! *** supermodule file not found! ");
177 return -1;
178 }
179 if (LoadSuperModuleFile(st2)) return -1;
180 }
181
483fd7d8 182 if (strstr(st,"SET_PARSIG_TRA")) {
183 sscanf(st,"%s %f",tmp,&f1);
184 fParSigTranslations=f1;
185 }
186
187 if (strstr(st,"SET_PARSIG_ROT")) {
188 sscanf(st,"%s %f",tmp,&f1);
189 fParSigRotations=f1;
190 }
191
192 if (strstr(st,"SET_NSTDDEV")) {
193 sscanf(st,"%s %d",tmp,&idx);
194 fNStdDev=idx;
195 }
196
197 if (strstr(st,"SET_RESCUT_INIT")) {
198 sscanf(st,"%s %f",tmp,&f1);
199 fResCutInitial=f1;
200 }
201
202 if (strstr(st,"SET_RESCUT_OTHER")) {
203 sscanf(st,"%s %f",tmp,&f1);
204 fResCut=f1;
205 }
206
207 if (strstr(st,"SET_STARTFAC")) {
208 sscanf(st,"%s %f",tmp,&f1);
209 fStartFac=f1;
210 }
211
212 if (strstr(st,"SET_LOCAL_SHIFTS")) { // only enabled mode...
213 fUseLocalShifts = kTRUE;
214 }
215
0856765c 216 if (strstr(st,"MODULE_INDEX")) { // works only for sensitive modules
483fd7d8 217 sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips);
218 voluid=GetModuleVolumeID(idx);
0856765c 219 if (!voluid || voluid>14300) return 1; // bad index
483fd7d8 220 fModuleIndex[nmod]=idx;
221 fModuleVolumeID[nmod]=voluid;
222 fFreeParam[nmod][0]=itx;
223 fFreeParam[nmod][1]=ity;
224 fFreeParam[nmod][2]=itz;
225 fFreeParam[nmod][3]=iph;
226 fFreeParam[nmod][4]=ith;
227 fFreeParam[nmod][5]=ips;
0856765c 228 fMilleModule[nmod] = new AliITSAlignMilleModule(voluid);
483fd7d8 229 nmod++;
230 }
231
232 if (strstr(st,"MODULE_VOLUID")) {
0856765c 233 sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips);
234 voluid=UShort_t(idx);
235 if (voluid>14335 && fUseSuperModules) { // custom supermodule
236 int ism=-1;
237 for (int j=0; j<fNSuperModules; j++) {
238 if (voluid==fSuperModule[j]->GetVolumeID()) ism=j;
239 }
240 if (ism<0) return -1; // bad volid
241 fModuleIndex[nmod]=fSuperModule[ism]->GetIndex();
242 fModuleVolumeID[nmod]=voluid;
243 fFreeParam[nmod][0]=itx;
244 fFreeParam[nmod][1]=ity;
245 fFreeParam[nmod][2]=itz;
246 fFreeParam[nmod][3]=iph;
247 fFreeParam[nmod][4]=ith;
248 fFreeParam[nmod][5]=ips;
249 fMilleModule[nmod] = new AliITSAlignMilleModule(*fSuperModule[ism]);
250 nmod++;
251 }
252 else { // sensitive volume
253 idx=GetModuleIndex(voluid);
254 if (idx<0 || idx>2197) return 1; // bad index
255 fModuleIndex[nmod]=idx;
256 fModuleVolumeID[nmod]=voluid;
257 fFreeParam[nmod][0]=itx;
258 fFreeParam[nmod][1]=ity;
259 fFreeParam[nmod][2]=itz;
260 fFreeParam[nmod][3]=iph;
261 fFreeParam[nmod][4]=ith;
262 fFreeParam[nmod][5]=ips;
263 fMilleModule[nmod] = new AliITSAlignMilleModule(voluid);
264 nmod++;
265 }
483fd7d8 266 }
0856765c 267 //----------
483fd7d8 268
269 } // end while
270
271 fNModules = nmod;
272 fNGlobal = fNModules*fgNParCh;
273
274 fclose(pfc);
275 return 0;
276}
277
278Int_t AliITSAlignMille::GetModuleIndex(const Char_t *symname) {
279 /// index from symname
280 if (!symname) return -1;
281 for (Int_t i=0; i<2198; i++) {
282 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
283 }
284 return -1;
285}
286
287Int_t AliITSAlignMille::GetModuleIndex(UShort_t voluid) {
288 /// index from volume ID
289 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
290 if (lay<1|| lay>6) return -1;
291 Int_t idx=Int_t(voluid)-2048*lay;
292 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
293 for (Int_t ilay=1; ilay<lay; ilay++)
294 idx += AliGeomManager::LayerSize(ilay);
295 return idx;
296}
297
298UShort_t AliITSAlignMille::GetModuleVolumeID(const Char_t *symname) {
299 /// volume ID from symname
0856765c 300 /// works for sensitive volumes only
483fd7d8 301 if (!symname) return 0;
302
303 for (UShort_t voluid=2000; voluid<13300; voluid++) {
304 Int_t modId;
305 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
306 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
307 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
308 }
309 }
310
311 return 0;
312}
313
314UShort_t AliITSAlignMille::GetModuleVolumeID(Int_t index) {
315 /// volume ID from index
0856765c 316 if (index<0) return 0;
317 if (index<2198)
318 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
319 else {
320 for (int i=0; i<fNSuperModules; i++) {
321 if (fSuperModule[i]->GetIndex()==index) return fSuperModule[i]->GetVolumeID();
322 }
323 }
324 return 0;
483fd7d8 325}
326
327void AliITSAlignMille::InitGeometry() {
328 /// initialize geometry
329 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
330 fGeoManager = AliGeomManager::GetGeometry();
331 if (!fGeoManager) {
332 AliInfo("Couldn't initialize geometry");
333 return;
334 }
335 // temporary align object, just use the rotation...
0856765c 336 //fTempAlignObj=new AliAlignObjParams(AliITSgeomTGeo::GetSymName(7),2055,0,0,0,0,0,0,kFALSE);
337 fTempAlignObj=new AliAlignObjParams;
483fd7d8 338}
339
340void AliITSAlignMille::Init(Int_t nGlobal, /* number of global paramers */
341 Int_t nLocal, /* number of local parameters */
342 Int_t nStdDev /* std dev cut */ )
343{
344 /// Initialization of AliMillepede. Fix parameters, define constraints ...
345 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
346 fIsMilleInit = kTRUE;
347
348 /// Fix non free parameters
349 for (Int_t i=0; i<fNModules; i++) {
350 for (Int_t j=0; j<ITSMILLE_NPARCH; j++) {
351 if (!fFreeParam[i][j]) FixParameter(i*ITSMILLE_NPARCH+j,0.0);
352 else {
353 // pepopepo: da sistemare il settaggio delle sigma...
354 Double_t parsig=0;
355 if (j<3) parsig=fParSigTranslations; // translations (0.0100 cm)
356 else parsig=fParSigRotations; // rotations (1/10 deg)
357 FixParameter(i*ITSMILLE_NPARCH+j,parsig);
358 }
359 }
360 }
361
362
363// // Set iterations
364 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
365}
366
367
368void AliITSAlignMille::AddConstraint(Double_t *par, Double_t value) {
369 /// Constrain equation defined by par to value
370 if (!fIsMilleInit) {
371 AliInfo("Millepede has not been initialized!");
372 return;
373 }
374 fMillepede->SetGlobalConstraint(par, value);
375 AliInfo("Adding constraint");
376}
377
378void AliITSAlignMille::InitGlobalParameters(Double_t *par) {
379 /// Initialize global parameters with par array
380 if (!fIsMilleInit) {
381 AliInfo("Millepede has not been initialized!");
382 return;
383 }
384 fMillepede->SetGlobalParameters(par);
385 AliInfo("Init Global Parameters");
386}
387
388void AliITSAlignMille::FixParameter(Int_t iPar, Double_t value) {
389 /// Parameter iPar is encourage to vary in [-value;value].
390 /// If value == 0, parameter is fixed
391 if (!fIsMilleInit) {
392 AliInfo("Millepede has not been initialized!");
393 return;
394 }
395 fMillepede->SetParSigma(iPar, value);
396 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
397}
398
399void AliITSAlignMille::ResetLocalEquation()
400{
401 /// Reset the derivative vectors
402 for(int i=0; i<fNLocal; i++) {
403 fLocalDerivatives[i] = 0.0;
404 }
405 for(int i=0; i<fNGlobal; i++) {
406 fGlobalDerivatives[i] = 0.0;
407 }
408}
409
0856765c 410Int_t AliITSAlignMille::ApplyToGeometry() {
411 /// apply starting realignment to ideal geometry
412 if (!fGeoManager) return -1;
413 TFile *pref = new TFile(fPreAlignmentFileName.Data());
414 if (!pref->IsOpen()) return -2;
415 TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
416 if (!prea) return -3;
417 Int_t nprea=prea->GetEntriesFast();
418 AliInfo(Form("Array of input misalignments with %d entries",nprea));
419
420 for (int ix=0; ix<nprea; ix++) {
421 AliAlignObjParams *preo=(AliAlignObjParams*) prea->UncheckedAt(ix);
422 if (!preo->ApplyToGeometry()) return -4;
423 }
424 pref->Close();
425 delete pref;
426
427 fUsePreAlignment = kTRUE;
428 return 0;
429}
430
483fd7d8 431Int_t AliITSAlignMille::InitModuleParams() {
432 /// initialize geometry parameters for a given detector
433 /// for current cluster (fCluster)
434 /// fGlobalInitParam[] is set as:
435 /// [tx,ty,tz,psi,theta,phi]
436 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
437 /// *** At the moment: using Raffalele's angles definition ***
438 ///
439 /// Num of Dets: ITSMILLE_NDETELEM = fgNDetElem
440 /// Num of Par : ITSMILLE_NPARCH = fgNParCh
441 /// return 0 if success
442
443 if (!fGeoManager) {
444 AliInfo("ITS geometry not initialized!");
445 return -1;
446 }
447
0856765c 448 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
449
483fd7d8 450 // set the internal index (index in module list)
451 UShort_t voluid=fCluster.GetVolumeID();
452 Int_t k=fNModules-1;
0856765c 453 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--; // new
483fd7d8 454 if (k<0) return -3;
0856765c 455 fCurrentModuleInternalIndex=k; // the internal index of the SUPERMODULE
456
457 fCurrentModuleIndex=fMilleModule[k]->GetIndex(); // index of the SUPERMODULE
483fd7d8 458
459 // set the index
0856765c 460 Int_t index = GetModuleIndex(voluid);
483fd7d8 461 if (index<0) return -2;
0856765c 462 fCurrentSensVolIndex = index; // the index of the SENSITIVE VOLUME
483fd7d8 463
464 fModuleInitParam[0] = 0.0;
465 fModuleInitParam[1] = 0.0;
466 fModuleInitParam[2] = 0.0;
467 fModuleInitParam[3] = 0.0; // psi (X)
468 fModuleInitParam[4] = 0.0; // theta (Y)
469 fModuleInitParam[5] = 0.0; // phi (Z)
470
471 /// get global (corrected) matrix
472 // if (!AliITSgeomTGeo::GetOrigMatrix(index,*fCurrentModuleHMatrix)) return -3;
0856765c 473// Double_t rott[9];
474// if (!AliITSgeomTGeo::GetRotation(index,rott)) return -3;
475// fCurrentModuleHMatrix->SetRotation(rott);
476// Double_t oLoc[3]={0,0,0};
477// if (!AliITSgeomTGeo::LocalToGlobal(index,oLoc,fCurrentModuleTranslation)) return -4;
478// fCurrentModuleHMatrix->SetTranslation(fCurrentModuleTranslation);
479
480// new
481 fCurrentModuleHMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetMatrix();
482
483 for (int ii=0; ii<3; ii++)
484 fCurrentModuleTranslation[ii]=fCurrentModuleHMatrix->GetTranslation()[ii];
485
486 TGeoHMatrix *svOrigMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeOrigGlobalMatrix(voluid);
483fd7d8 487
488 /// get back local coordinates
489 fMeasGlo[0] = fCluster.GetX();
490 fMeasGlo[1] = fCluster.GetY();
491 fMeasGlo[2] = fCluster.GetZ();
0856765c 492 svOrigMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
493 //svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
494 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
495
496 TGeoHMatrix *svMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeMatrix(voluid);
497
498 // modify global coordinates according with pre-aligment
499 svMatrix->LocalToMaster(fMeasLoc,fMeasGlo);
500 fCluster.SetXYZ(fMeasGlo[0],fMeasGlo[1] ,fMeasGlo[2]);
501 AliDebug(2,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasGlo[0] ,fMeasGlo[1] ,fMeasGlo[2] ));
483fd7d8 502
0856765c 503 // mettere il new GetLocalSigma...
483fd7d8 504 // set stdev from cluster
505 TGeoHMatrix hcov;
506 Double_t hcovel[9];
507 hcovel[0]=double(fCluster.GetCov()[0]);
508 hcovel[1]=double(fCluster.GetCov()[1]);
509 hcovel[2]=double(fCluster.GetCov()[3]);
510 hcovel[3]=double(fCluster.GetCov()[1]);
511 hcovel[4]=double(fCluster.GetCov()[2]);
512 hcovel[5]=double(fCluster.GetCov()[4]);
513 hcovel[6]=double(fCluster.GetCov()[3]);
514 hcovel[7]=double(fCluster.GetCov()[4]);
515 hcovel[8]=double(fCluster.GetCov()[5]);
516 hcov.SetRotation(hcovel);
517 // now rotate in local system
0856765c 518 hcov.MultiplyLeft(&svMatrix->Inverse());
519 hcov.Multiply(svMatrix);
483fd7d8 520
521 // per i ruotati c'e' delle sigmaY che compaiono... prob
522 // e' un problema di troncamento
523 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
524 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4]));
525 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
526
0856765c 527 // set minimum value for SigmaLoc to 10 micron
528 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
529 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
530
483fd7d8 531 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%f fSigmaLocY=%f fSigmaLocZ=%f \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
532
533 return 0;
534}
535
536void AliITSAlignMille::SetCurrentModule(Int_t index) {
0856765c 537 /// set as current the SuperModule that contains the 'index' sens.vol.
538 if (index<0 || index>2197) {
539 AliInfo("index does not correspond to a sensitive volume!");
540 return;
541 }
483fd7d8 542 UShort_t voluid=GetModuleVolumeID(index);
0856765c 543 //Int_t k=IsDefined(voluid);
544 Int_t k=IsContained(voluid);
545 if (k>=0){
546 //if (voluid<14336)
547 fCluster.SetVolumeID(voluid);
548 //else {
549 //fCluster.SetVolumeID(fMilleModule[k]->GetSensitiveVolumeVolumeID()[0]);
550 //printf("current module is a supermodule: fCluster set to first sensitive volume of the supermodule\n");
551 //}
552 fCluster.SetXYZ(0,0,0);
553 InitModuleParams();
554 }
555 else
556 printf("module %d not defined\n",index);
557}
558
559void AliITSAlignMille::SetCurrentSensitiveModule(Int_t index) {
560 /// set as current the SuperModule that contains the 'index' sens.vol.
561 if (index<0 || index>2197) {
562 AliInfo("index does not correspond to a sensitive volume!");
563 return;
564 }
565 UShort_t voluid=AliITSAlignMilleModule::GetVolumeIDFromIndex(index);
566 Int_t k=IsDefined(voluid);
567 //printf("---> voluid=%d k=%d\n",voluid,k);
568 if (k>=0){
483fd7d8 569 fCluster.SetVolumeID(voluid);
570 fCluster.SetXYZ(0,0,0);
571 InitModuleParams();
572 }
0856765c 573 else
574 printf("module %d not defined\n",index);
483fd7d8 575}
576
0856765c 577void AliITSAlignMille::Print(Option_t*) const
578{
483fd7d8 579 ///
580 printf("*** AliMillepede for ITS ***\n");
0856765c 581 printf(" number of defined super modules: %d\n",fNModules);
582
483fd7d8 583 if (fGeoManager)
584 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
585 else
586 printf(" geometry not loaded\n");
0856765c 587
588 if (fUseSuperModules)
589 printf(" using custom supermodules ( %d defined )\n",fNSuperModules);
590 else
591 printf(" custom supermodules not used\n");
592
593 if (fUsePreAlignment)
594 printf(" using prealignment from %s \n",fPreAlignmentFileName.Data());
595 else
596 printf(" prealignment not used\n");
597
483fd7d8 598 if (fUseLocalShifts)
599 printf(" Alignment shifts will be computed in LOCAL RS\n");
600 else
601 printf(" Alignment shifts will be computed in GLOBAL RS\n");
0856765c 602
483fd7d8 603 printf(" Millepede configuration parameters:\n");
604 printf(" init parsig for translations : %.4f\n",fParSigTranslations);
605 printf(" init parsig for rotations : %.4f\n",fParSigRotations);
606 printf(" init value for chi2 cut : %.4f\n",fStartFac);
607 printf(" first iteration cut value : %.4f\n",fResCutInitial);
608 printf(" other iterations cut value : %.4f\n",fResCut);
609 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
610
611 printf("List of defined modules:\n");
0856765c 612 printf(" intidx\tindex\tvoluid\tname\n");
483fd7d8 613 for (int i=0; i<fNModules; i++)
0856765c 614 printf(" %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],fMilleModule[i]->GetName());
615
483fd7d8 616}
617
0856765c 618AliITSAlignMilleModule *AliITSAlignMille::GetMilleModule(UShort_t voluid)
619{
620 // return pointer to a define supermodule
621 // return NULL if error
622 Int_t i=IsDefined(voluid);
623 if (i<0) return NULL;
624 return fMilleModule[i];
625}
626
627AliITSAlignMilleModule *AliITSAlignMille::GetCurrentModule()
628{
629 if (fNModules) return fMilleModule[fCurrentModuleInternalIndex];
630 return NULL;
631}
632
633void AliITSAlignMille::PrintCurrentModuleInfo()
634{
483fd7d8 635 ///
0856765c 636 Int_t k=fCurrentModuleInternalIndex;
637 if (k<0 || k>=fNModules) return;
638 fMilleModule[k]->Print("");
483fd7d8 639}
640
641
642void AliITSAlignMille::InitTrackParams(int meth) {
643 /// initialize local parameters with different methods
644 /// for current track (fTrack)
645
0856765c 646 Int_t npts=0;
647 TF1 *f1=NULL;
648 TGraph *g=NULL;
649 Float_t sigmax[20],sigmay[20],sigmaz[20];
650 AliTrackPoint ap;
651 TGraphErrors *ge=NULL;
652
483fd7d8 653 switch (meth) {
654 case 1: // simple linear interpolation
655 // get local starting parameters (to be substituted by ESD track parms)
656 // local parms (fLocalInitParam[]) are:
657 // [0] = global x coord. of straight line intersection at y=0 plane
658 // [1] = global z coord. of straight line intersection at y=0 plane
659 // [2] = px/py
660 // [3] = pz/py
661
662 // test #1: linear fit in x(y) and z(y)
0856765c 663 npts = fTrack->GetNPoints();
483fd7d8 664
0856765c 665 f1=new TF1("f1","[0]+x*[1]",-50,50);
483fd7d8 666
0856765c 667 g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
483fd7d8 668 g->Fit(f1,"RNQ");
669 fLocalInitParam[0] = f1->GetParameter(0);
670 fLocalInitParam[2] = f1->GetParameter(1);
671 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
672 delete g; g=NULL;
673
674 g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ());
675 g->Fit(f1,"RNQ");
676 fLocalInitParam[1] = f1->GetParameter(0);
677 fLocalInitParam[3] = f1->GetParameter(1);
678 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
679 delete g;
680 delete f1;
681
682 break;
0856765c 683
684 case 2: // simple linear interpolation weighted using sigmas
685 // get local starting parameters (to be substituted by ESD track parms)
686 // local parms (fLocalInitParam[]) are:
687 // [0] = global x coord. of straight line intersection at y=0 plane
688 // [1] = global z coord. of straight line intersection at y=0 plane
689 // [2] = px/py
690 // [3] = pz/py
691
692 // test #1: linear fit in x(y) and z(y)
693 npts = fTrack->GetNPoints();
694 for (Int_t isig=0; isig<npts; isig++) {
695 fTrack->GetPoint(ap,isig);
696 sigmax[isig]=ap.GetCov()[0];
697 if (sigmax[isig]<1.0e-07) sigmax[isig]=1.0e-07; // minimum sigma=3 mu
698 sigmax[isig]=TMath::Sqrt(sigmax[isig]);
699
700 sigmay[isig]=ap.GetCov()[2];
701 if (sigmay[isig]<1.0e-07) sigmay[isig]=1.0e-07; // minimum sigma=3 mu
702 sigmay[isig]=TMath::Sqrt(sigmay[isig]);
703
704 sigmaz[isig]=ap.GetCov()[5];
705 if (sigmaz[isig]<1.0e-07) sigmaz[isig]=1.0e-07; // minimum sigma=3 mu
706 sigmaz[isig]=TMath::Sqrt(sigmaz[isig]);
707 }
708
709 f1=new TF1("f1","[0]+x*[1]",-50,50);
710
711 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetX(),sigmay,sigmax);
712 ge->Fit(f1,"RNQ");
713 fLocalInitParam[0] = f1->GetParameter(0);
714 fLocalInitParam[2] = f1->GetParameter(1);
715 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
716 delete ge; ge=NULL;
717
718 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetZ(),sigmay,sigmaz);
719 ge->Fit(f1,"RNQ");
720 fLocalInitParam[1] = f1->GetParameter(0);
721 fLocalInitParam[3] = f1->GetParameter(1);
722 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
723 delete ge;
724 delete f1;
725
726 break;
727
483fd7d8 728 }
0856765c 729}
483fd7d8 730
0856765c 731Int_t AliITSAlignMille::IsDefined(UShort_t voluid) const
732{
733 // checks if supermodule 'voluid' is defined and return the internal index
734 // return -1 if error
735 Int_t k=fNModules-1;
736 while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
737 if (k<0) return -1;
738 return k;
739}
740
741Int_t AliITSAlignMille::IsContained(UShort_t voluid) const
742{
743 // checks if the sensitive module 'voluid' is contained inside a supermodule and return the internal index of the last identified supermodule
744 // return -1 if error
745 if (AliITSAlignMilleModule::GetIndexFromVolumeID(voluid)<0) return -1;
746 Int_t k=fNModules-1;
747 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
748 if (k<0) return -1;
749 return k;
483fd7d8 750}
0856765c 751
483fd7d8 752Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
753{
0856765c 754 /// check if a sensitive volume is contained inside one of the defined supermodules
483fd7d8 755 Int_t k=fNModules-1;
0856765c 756 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
483fd7d8 757 if (k>=0) return kTRUE;
758 return kFALSE;
759}
760
761Int_t AliITSAlignMille::CheckCurrentTrack() {
762 /// checks if AliTrackPoints belongs to defined modules
763 /// return number of good poins
764 /// return 0 if not enough points
765
766 Int_t npts = fTrack->GetNPoints();
767 Int_t ngoodpts=0;
768 // debug points
769 for (int j=0; j<npts; j++) {
770 UShort_t voluid = fTrack->GetVolumeID()[j];
771 if (CheckVolumeID(voluid)) {
772 ngoodpts++;
773 }
774 }
775 // pepo da controllare...
776 if (ngoodpts<fMinNPtsPerTrack) return 0;
777
778 return ngoodpts;
779}
780
781Int_t AliITSAlignMille::ProcessTrack(AliTrackPointArray *track) {
782 /// Process track; Loop over hits and set local equations
783 /// here 'track' is a AliTrackPointArray
784 /// return 0 if success;
785
786 if (!fIsMilleInit) {
787 AliInfo("Millepede has not been initialized!");
788 return -1;
789 }
790
791 Int_t npts = track->GetNPoints();
792 AliDebug(2,Form("\n*** Processing track with %d points ***\n",npts));
793 fTrack = track;
794
795 // checks if there are enough good points
796 if (!CheckCurrentTrack()) {
797 AliInfo("Track with not enough good points - will not be used...");
798 return -1;
799 }
800
801 // set local starting parameters (to be substituted by ESD track parms)
802 // local parms (fLocalInitParam[]) are:
803 // [0] = global x coord. of straight line intersection at y=0 plane
804 // [1] = global z coord. of straight line intersection at y=0 plane
805 // [2] = px/py
806 // [3] = pz/py
0856765c 807 InitTrackParams(fInitTrackParamsMeth);
483fd7d8 808
809 for (Int_t ipt=0; ipt<npts; ipt++) {
810 fTrack->GetPoint(fCluster,ipt);
811 if (!CheckVolumeID(fCluster.GetVolumeID())) continue;
0856765c 812 AliDebug(2,Form(" Original Point = ( %f , %f , %f ) volid=%d\n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ(),fCluster.GetVolumeID()));
483fd7d8 813
814 // set geometry parameters for the the current module
815 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
816 if (InitModuleParams()) continue;
817
818 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
0856765c 819 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
483fd7d8 820
821 if (SetLocalEquations()) return -1;
822
823 } // end loop over points
824
825 return 0;
826}
827
828Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) {
829 /// calculate track intersection point in local coordinates
830 /// according with a given set of parameters (local(4) and global(6))
831 /// and fill fPintLoc/Glo
832 /// local are: pgx0, pgz0, ugx0, ugz0
833 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
834 /// return 0 if success
835
836 AliDebug(3,Form("lpar = %g %g %g %g \ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
837
838 // vector of initial straight line direction in glob. coord
839 // ATTENZIONE: forse va rinormalizzato tutto...
840 Double_t v0g[3];
841 //Double_t
842 v0g[0]=lpar[2];
843 v0g[1]=1.0;
844 v0g[2]=lpar[3];
845
846 // intercept in yg=0 plane in glob coord
847 Double_t p0g[3];
848 p0g[0]=lpar[0];
849 p0g[1]=0.0;
850 p0g[2]=lpar[1];
851
852
853 // prepare the TGeoHMatrix
0856765c 854// Double_t tr[3],ang[3];
855// //Double_t rad2deg=180./TMath::Pi();
856// if (fUseLocalShifts) { // just Delta matrix
857// tr[0]=gpar[0];
858// tr[1]=gpar[1];
859// tr[2]=gpar[2];
860// ang[0]=gpar[3]; // psi (X)
861// ang[1]=gpar[4]; // theta (Y)
862// ang[2]=gpar[5]; // phi (Z)
863// }
864// else { // total matrix with shifted parameter
865// AliInfo("global shifts not implemented yet!");
866// return -1;
867// }
868
869// //printf("fTempRot = 0x%x - ang = %g %g %g \n",fTempRot,gpar[5]*rad2deg,gpar[3]*rad2deg,gpar[4]*rad2deg);
870
871// fTempAlignObj->SetRotation(ang[0],ang[1],ang[2]);
872// AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2]));
873// TGeoHMatrix hm;
874// fTempAlignObj->GetMatrix(hm);
875// fTempHMat->SetRotation(hm.GetRotationMatrix());
876// fTempHMat->SetTranslation(tr);
483fd7d8 877
0856765c 878// // in this case the gpar[] array contains only shifts
879// // and fInitModuleParam[] are set to 0
880// // fCurrentModuleHMatrix is then modified as fCurrentHM*fTempHM
881// if (fUseLocalShifts)
882// fTempHMat->MultiplyLeft(fCurrentModuleHMatrix);
883 TGeoHMatrix *fTempHMat = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar);
884 if (!fTempHMat) return -1;
483fd7d8 885
886 // same in local coord.
887 Double_t p0l[3],v0l[3];
888 fTempHMat->MasterToLocalVect(v0g,v0l);
889 fTempHMat->MasterToLocal(p0g,p0l);
890
891 if (TMath::Abs(v0l[1])<1e-15) {
892 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
893 return -1;
894 }
895
896 // local intersection point
897 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
898 fPintLoc[1] = 0;
899 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
900
901 // global intersection point
902 fTempHMat->LocalToMaster(fPintLoc,fPintGlo);
903 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]));
904
905 return 0;
906}
907
908Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) {
909 /// calculate numerically (ROOT's style) the derivatives for
910 /// local X intersection and local Z intersection
911 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0
912 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
913 /// return 0 if success
914
915 // copy initial parameters
916 Double_t lpar[ITSMILLE_NLOCAL];
917 Double_t gpar[ITSMILLE_NPARCH];
918 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) lpar[i]=fLocalInitParam[i];
919 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) gpar[i]=fModuleInitParam[i];
920
921 // pepopepo
922 // trial with fixed dpar...
923 Double_t dpar=0.0;
924 if (islpar) {
925 //dpar=fLocalInitParam[paridx]*0.001;
926 // set minimum dpar
927 if (paridx<2) dpar=1.0e-4; // translations
928 else dpar=1.0e-6; // direction
929 }
930 else {
931 //dpar=fModuleInitParam[paridx]*0.001;
932 if (paridx<3) dpar=1.0e-4; // translations
933 else dpar=1.0e-2; // angles
934 }
935 AliDebug(3,Form("\n+++ automatic dpar=%g\n",dpar));
936 if (fDeltaPar) dpar=fDeltaPar;
937 AliDebug(3,Form("+++ using dpar=%g\n\n",dpar));
938
939 // calculate derivative ROOT's like:
940 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
941 Double_t pintl1[3]; // f(x-h)
942 Double_t pintl2[3]; // f(x-h/2)
943 Double_t pintl3[3]; // f(x+h/2)
944 Double_t pintl4[3]; // f(x+h)
945
946 // first values
947 if (islpar) lpar[paridx] -= dpar;
948 else gpar[paridx] -= dpar;
949 if (CalcIntersectionPoint(lpar, gpar)) return -2;
950 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
951
952 // second values
953 if (islpar) lpar[paridx] += dpar/2;
954 else gpar[paridx] += dpar/2;
955 if (CalcIntersectionPoint(lpar, gpar)) return -2;
956 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
957
958 // third values
959 if (islpar) lpar[paridx] += dpar;
960 else gpar[paridx] += dpar;
961 if (CalcIntersectionPoint(lpar, gpar)) return -2;
962 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
963
964 // fourth values
965 if (islpar) lpar[paridx] += dpar/2;
966 else gpar[paridx] += dpar/2;
967 if (CalcIntersectionPoint(lpar, gpar)) return -2;
968 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
969
970 Double_t h2 = 1./(2.*dpar);
971 Double_t d0 = pintl4[0]-pintl1[0];
972 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
973 fDerivativeXLoc = h2*(4*d2 - d0)/3.;
974 if (TMath::Abs(fDerivativeXLoc) < 1.0e-9) fDerivativeXLoc=0.0;
975
976 d0 = pintl4[2]-pintl1[2];
977 d2 = 2.*(pintl3[2]-pintl2[2]);
978 fDerivativeZLoc = h2*(4*d2 - d0)/3.;
979 if (TMath::Abs(fDerivativeZLoc) < 1.0e-9) fDerivativeZLoc=0.0;
980
981 AliDebug(3,Form("\n+++ derivatives +++ \n"));
982 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",fDerivativeXLoc));
983 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",fDerivativeZLoc));
984
985 return 0;
986}
987
988Int_t AliITSAlignMille::SetLocalEquations() {
989 /// Define local equation for current track and cluster in x coor.
990 /// return 0 if success
991
992 // store first interaction point
993 CalcIntersectionPoint(fLocalInitParam, fModuleInitParam);
994 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
0856765c 995 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
483fd7d8 996
997 // calculate local derivatives numerically
998 Double_t dXdL[ITSMILLE_NLOCAL],dZdL[ITSMILLE_NLOCAL];
999 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) {
1000 if (CalcDerivatives(i,kTRUE)) return -1;
1001 dXdL[i]=fDerivativeXLoc;
1002 dZdL[i]=fDerivativeZLoc;
1003 }
1004
1005 Double_t dXdG[ITSMILLE_NPARCH],dZdG[ITSMILLE_NPARCH];
1006 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
1007 if (CalcDerivatives(i,kFALSE)) return -1;
1008 dXdG[i]=fDerivativeXLoc;
1009 dZdG[i]=fDerivativeZLoc;
1010 }
1011
1012 AliDebug(2,Form("\n***************\n"));
1013 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
1014 AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
1015 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1016 AliDebug(2,Form("Global parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdG[i],dZdG[i]));
1017 AliDebug(2,Form("\n***************\n"));
1018
1019
1020 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
1021 // set equation for Xloc coordinate
1022 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
1023 SetLocalDerivative(i,dXdL[i]);
1024 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1025 SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dXdG[i]);
1026 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]);
1027
1028
1029 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
1030 // set equation for Zloc coordinate
1031 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
1032 SetLocalDerivative(i,dZdL[i]);
1033 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1034 SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dZdG[i]);
1035 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]);
1036
1037 return 0;
1038}
1039
1040
1041void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
1042 /// Call local fit for this track
1043 if (!fIsMilleInit) {
1044 AliInfo("Millepede has not been initialized!");
1045 return;
1046 }
1047 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
1048 AliDebug(2,Form("iRes = %d",iRes));
1049 if (iRes && !lSingleFit) {
1050 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
1051 }
1052}
1053
1054void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
1055 /// Call global fit; Global parameters are stored in parameters
1056 if (!fIsMilleInit) {
1057 AliInfo("Millepede has not been initialized!");
1058 return;
1059 }
1060 fMillepede->GlobalFit(parameters,errors,pulls);
1061 AliInfo("Done fitting global parameters!");
1062}
1063
1064Double_t AliITSAlignMille::GetParError(Int_t iPar) {
1065 /// Get error of parameter iPar
1066 if (!fIsMilleInit) {
1067 AliInfo("Millepede has not been initialized!");
1068 return 0;
1069 }
1070 Double_t lErr = fMillepede->GetParError(iPar);
1071 return lErr;
1072}
1073
1074void AliITSAlignMille::PrintGlobalParameters() {
1075 /// Print global parameters
1076 if (!fIsMilleInit) {
1077 AliInfo("Millepede has not been initialized!");
1078 return;
1079 }
1080 fMillepede->PrintGlobalParameters();
1081}
1082
1083// //_________________________________________________________________________
0856765c 1084Int_t AliITSAlignMille::LoadSuperModuleFile(const Char_t *sfile)
1085{
1086 // load definitions of supermodules from a root file
1087 // return 0 if success
1088
1089 TFile *smf=TFile::Open(sfile);
1090 if (!smf->IsOpen()) {
1091 AliInfo(Form("Cannot open supermodule file %s",sfile));
1092 return -1;
1093 }
1094
1095 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
1096 if (!sma) {
1097 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
1098 return -2;
1099 }
1100 Int_t nsma=sma->GetEntriesFast();
1101 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
1102
1103 Char_t st[250];
1104 char symname[150];
1105 UShort_t volid;
1106 TGeoHMatrix m;
1107
1108 for (Int_t i=0; i<nsma; i++) {
1109 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
1110 volid=a->GetVolUID();
1111 strcpy(st,a->GetSymName());
1112 a->GetMatrix(m);
1113
1114 sscanf(st,"%s",symname);
1115 // decode module list
1116 char *stp=strstr(st,"ModuleList:");
1117 if (!stp) return -3;
1118 stp += 11;
1119 int idx[2200];
1120 char spp[200]; int jp=0;
1121 char cl[20];
1122 strcpy(st,stp);
1123 int l=strlen(st);
1124 int j=0;
1125 int n=0;
1126
1127 while (j<=l) {
1128 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
1129 spp[jp]=0;
1130 jp=0;
1131 if (strlen(spp)) {
1132 int k=strcspn(spp,"-");
1133 if (k<int(strlen(spp))) { // c'e' il -
1134 strcpy(cl,&(spp[k+1]));
1135 spp[k]=0;
1136 int ifrom=atoi(spp); int ito=atoi(cl);
1137 for (int b=ifrom; b<=ito; b++) {
1138 idx[n]=b;
1139 n++;
1140 }
1141 }
1142 else { // numerillo singolo
1143 idx[n]=atoi(spp);
1144 n++;
1145 }
1146 }
1147 }
1148 else {
1149 spp[jp]=st[j];
1150 jp++;
1151 }
1152 j++;
1153 }
1154 UShort_t volidsv[2198];
1155 for (j=0;j<n;j++) {
1156 volidsv[j]=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx[j]);
1157 if (!volidsv[j]) {
1158 AliInfo(Form("Index %d not valid (range 0->2197)",idx[j]));
1159 return -5;
1160 }
1161 }
1162 Int_t smindex=int(2198+volid-14336); // virtual index
1163 fSuperModule[fNSuperModules]=new AliITSAlignMilleModule(smindex,volid,symname,&m,n,volidsv);
1164
1165 //-------------
1166 fNSuperModules++;
1167 }
1168
1169 smf->Close();
1170
1171 fUseSuperModules=1;
1172 return 0;
1173}
483fd7d8 1174