]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSAlignMille.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Simple implementation of the VEventPool which allows to loop several times over the...
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille.cxx
... / ...
CommitLineData
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>
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
38#include "AliITSAlignMilleModule.h"
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),
70 fTempAlignObj(NULL),
71 fDerivativeXLoc(0),
72 fDerivativeZLoc(0),
73 fDeltaPar(0),
74 fMinNPtsPerTrack(3),
75 fInitTrackParamsMeth(1),
76 fGeometryFileName("geometry.root"),
77 fPreAlignmentFileName(""),
78 fGeoManager(0),
79 fCurrentModuleIndex(0),
80 fCurrentModuleInternalIndex(0),
81 fCurrentSensVolIndex(0),
82 fNModules(0),
83 fUseLocalShifts(kTRUE),
84 fUseSuperModules(kFALSE),
85 fUsePreAlignment(kFALSE),
86 fNSuperModules(0),
87 fCurrentModuleHMatrix(NULL)
88{
89 /// main constructor that takes input from configuration file
90
91 fMillepede = new AliMillepede();
92 fGlobalDerivatives = new Double_t[fNGlobal];
93 //fTempHMat = new TGeoHMatrix;
94 //fCurrentModuleHMatrix = new TGeoHMatrix;
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;
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];
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
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
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
216 if (strstr(st,"MODULE_INDEX")) { // works only for sensitive modules
217 sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips);
218 voluid=GetModuleVolumeID(idx);
219 if (!voluid || voluid>14300) return 1; // bad index
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;
228 fMilleModule[nmod] = new AliITSAlignMilleModule(voluid);
229 nmod++;
230 }
231
232 if (strstr(st,"MODULE_VOLUID")) {
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 }
266 }
267 //----------
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
300 /// works for sensitive volumes only
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
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;
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...
336 //fTempAlignObj=new AliAlignObjParams(AliITSgeomTGeo::GetSymName(7),2055,0,0,0,0,0,0,kFALSE);
337 fTempAlignObj=new AliAlignObjParams;
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
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
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
448 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
449
450 // set the internal index (index in module list)
451 UShort_t voluid=fCluster.GetVolumeID();
452 Int_t k=fNModules-1;
453 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--; // new
454 if (k<0) return -3;
455 fCurrentModuleInternalIndex=k; // the internal index of the SUPERMODULE
456
457 fCurrentModuleIndex=fMilleModule[k]->GetIndex(); // index of the SUPERMODULE
458
459 // set the index
460 Int_t index = GetModuleIndex(voluid);
461 if (index<0) return -2;
462 fCurrentSensVolIndex = index; // the index of the SENSITIVE VOLUME
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;
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);
487
488 /// get back local coordinates
489 fMeasGlo[0] = fCluster.GetX();
490 fMeasGlo[1] = fCluster.GetY();
491 fMeasGlo[2] = fCluster.GetZ();
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] ));
502
503 // mettere il new GetLocalSigma...
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
518 hcov.MultiplyLeft(&svMatrix->Inverse());
519 hcov.Multiply(svMatrix);
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
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
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) {
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 }
542 UShort_t voluid=GetModuleVolumeID(index);
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){
569 fCluster.SetVolumeID(voluid);
570 fCluster.SetXYZ(0,0,0);
571 InitModuleParams();
572 }
573 else
574 printf("module %d not defined\n",index);
575}
576
577void AliITSAlignMille::Print(Option_t*) const
578{
579 ///
580 printf("*** AliMillepede for ITS ***\n");
581 printf(" number of defined super modules: %d\n",fNModules);
582
583 if (fGeoManager)
584 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
585 else
586 printf(" geometry not loaded\n");
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
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");
602
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");
612 printf(" intidx\tindex\tvoluid\tname\n");
613 for (int i=0; i<fNModules; i++)
614 printf(" %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],fMilleModule[i]->GetName());
615
616}
617
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{
635 ///
636 Int_t k=fCurrentModuleInternalIndex;
637 if (k<0 || k>=fNModules) return;
638 fMilleModule[k]->Print("");
639}
640
641
642void AliITSAlignMille::InitTrackParams(int meth) {
643 /// initialize local parameters with different methods
644 /// for current track (fTrack)
645
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
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)
663 npts = fTrack->GetNPoints();
664
665 f1=new TF1("f1","[0]+x*[1]",-50,50);
666
667 g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
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;
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
728 }
729}
730
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;
750}
751
752Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
753{
754 /// check if a sensitive volume is contained inside one of the defined supermodules
755 Int_t k=fNModules-1;
756 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
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
807 InitTrackParams(fInitTrackParamsMeth);
808
809 for (Int_t ipt=0; ipt<npts; ipt++) {
810 fTrack->GetPoint(fCluster,ipt);
811 if (!CheckVolumeID(fCluster.GetVolumeID())) continue;
812 AliDebug(2,Form(" Original Point = ( %f , %f , %f ) volid=%d\n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ(),fCluster.GetVolumeID()));
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]) ));
819 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
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
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);
877
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;
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];
995 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
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// //_________________________________________________________________________
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}
1174