]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSAlignMille.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Possibility of using tracks in magnetic field (with AliTrackFitterRiemann) - M. Lunardon
[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
d3603b4d 20/// Alignment class for the ALICE ITS detector
483fd7d8 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"
d3603b4d 46#include "TSystem.h"
47#include "AliTrackFitterRieman.h"
483fd7d8 48
49/// \cond CLASSIMP
50ClassImp(AliITSAlignMille)
51/// \endcond
52
53Int_t AliITSAlignMille::fgNDetElem = ITSMILLE_NDETELEM;
54Int_t AliITSAlignMille::fgNParCh = ITSMILLE_NPARCH;
55
56AliITSAlignMille::AliITSAlignMille(const Char_t *configFilename, Bool_t initmille)
57 : TObject(),
58 fMillepede(0),
59 fStartFac(16.),
60 fResCutInitial(100.),
61 fResCut(100.),
62 fNGlobal(ITSMILLE_NDETELEM*ITSMILLE_NPARCH),
d3603b4d 63 fNLocal(4),
483fd7d8 64 fNStdDev(ITSMILLE_NSTDEV),
65 fIsMilleInit(kFALSE),
66 fParSigTranslations(0.0100),
67 fParSigRotations(0.1),
68 fTrack(NULL),
69 fCluster(),
70 fGlobalDerivatives(NULL),
d3603b4d 71 fSigmaXfactor(1.0),
72 fSigmaZfactor(1.0),
483fd7d8 73 fTempAlignObj(NULL),
74 fDerivativeXLoc(0),
75 fDerivativeZLoc(0),
483fd7d8 76 fMinNPtsPerTrack(3),
0856765c 77 fInitTrackParamsMeth(1),
d3603b4d 78 fProcessedPoints(NULL),
79 fTotBadLocEqPoints(0),
80 fRieman(NULL),
81 fRequirePoints(kFALSE),
82 fTempExcludedModule(-1),
483fd7d8 83 fGeometryFileName("geometry.root"),
0856765c 84 fPreAlignmentFileName(""),
483fd7d8 85 fGeoManager(0),
86 fCurrentModuleIndex(0),
87 fCurrentModuleInternalIndex(0),
0856765c 88 fCurrentSensVolIndex(0),
483fd7d8 89 fNModules(0),
90 fUseLocalShifts(kTRUE),
0856765c 91 fUseSuperModules(kFALSE),
92 fUsePreAlignment(kFALSE),
d3603b4d 93 fUseSortedTracks(kTRUE),
94 fBOn(kFALSE),
95 fBField(0.0),
0856765c 96 fNSuperModules(0),
d3603b4d 97 fCurrentModuleHMatrix(NULL),
98 fIsConfigured(kFALSE),
99 fBug(0)
483fd7d8 100{
101 /// main constructor that takes input from configuration file
102
103 fMillepede = new AliMillepede();
104 fGlobalDerivatives = new Double_t[fNGlobal];
d3603b4d 105
106 for (Int_t i=0; i<ITSMILLE_NDETELEM*2; i++) {
107 fPreAlignQF[i]=-1;
108 fSensVolSigmaXfactor[i]=1.0;
109 fSensVolSigmaZfactor[i]=1.0;
110 }
111
112 for (Int_t i=0; i<6; i++) {
113 fNReqLayUp[i]=0;
114 fNReqLayDown[i]=0;
115 fNReqLay[i]=0;
116 }
117 for (Int_t i=0; i<3; i++) {
118 fNReqDetUp[i]=0;
119 fNReqDetDown[i]=0;
120 fNReqDet[i]=0;
121 }
122
483fd7d8 123 Int_t lc=LoadConfig(configFilename);
124 if (lc) {
125 AliInfo(Form("Error %d loading configuration from %s",lc,configFilename));
126 }
127 else {
d3603b4d 128 fIsConfigured=kTRUE;
483fd7d8 129 if (initmille && fNGlobal<20000) {
130 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
131 Init(fNGlobal, fNLocal, fNStdDev);
132 ResetLocalEquation();
133 AliInfo("Parameters initialized to zero");
134 }
135 else {
136 AliInfo("Millepede has not been initialized ... ");
137 }
138 }
139
d3603b4d 140 if (fNModules) {
141 fProcessedPoints=new Int_t[fNModules];
142 for (Int_t ipp=0; ipp<fNModules; ipp++) fProcessedPoints[ipp]=0;
143 }
483fd7d8 144}
145
146AliITSAlignMille::~AliITSAlignMille() {
147 /// Destructor
148 if (fMillepede) delete fMillepede;
149 delete [] fGlobalDerivatives;
0856765c 150 for (int i=0; i<fNModules; i++) delete fMilleModule[i];
151 for (int i=0; i<fNSuperModules; i++) delete fSuperModule[i];
d3603b4d 152 if (fNModules) delete [] fProcessedPoints;
153 if (fRieman) delete fRieman;
483fd7d8 154}
155
156///////////////////////////////////////////////////////////////////////
483fd7d8 157Int_t AliITSAlignMille::LoadConfig(const Char_t *cfile) {
158 /// return 0 if success
159 /// 1 if error in module index or voluid
160
161 FILE *pfc=fopen(cfile,"r");
162 if (!pfc) return -1;
163
164 Char_t st[200],st2[200];
165 Char_t tmp[100];
166 Int_t idx,itx,ity,itz,ith,ips,iph;
d3603b4d 167 Float_t f1,f2;
483fd7d8 168 UShort_t voluid;
169 Int_t nmod=0;
170
171 while (fgets(st,200,pfc)) {
172
173 // skip comments
174 for (int i=0; i<int(strlen(st)); i++) {
175 if (st[i]=='#') st[i]=0;
176 }
177
178 if (strstr(st,"GEOMETRY_FILE")) {
179 sscanf(st,"%s %s",tmp,st2);
180 if (gSystem->AccessPathName(st2)) {
181 AliInfo("*** WARNING! *** geometry file not found! ");
182 return -1;
183 }
184 fGeometryFileName=st2;
185 InitGeometry();
186 }
187
0856765c 188 if (strstr(st,"PREALIGNMENT_FILE")) {
189 sscanf(st,"%s %s",tmp,st2);
190 if (gSystem->AccessPathName(st2)) {
191 AliInfo("*** WARNING! *** prealignment file not found! ");
192 return -1;
193 }
194 fPreAlignmentFileName=st2;
195 itx=ApplyToGeometry();
196 if (itx) {
197 AliInfo(Form("*** WARNING! *** error %d reading prealignment file! ",itx));
198 return -6;
199 }
200 }
201
202 if (strstr(st,"SUPERMODULE_FILE")) {
203 sscanf(st,"%s %s",tmp,st2);
204 if (gSystem->AccessPathName(st2)) {
205 AliInfo("*** WARNING! *** supermodule file not found! ");
206 return -1;
207 }
208 if (LoadSuperModuleFile(st2)) return -1;
209 }
210
d3603b4d 211 if (strstr(st,"SET_B_FIELD")) {
212 sscanf(st,"%s %f",tmp,&f1);
213 if (f1>0) {
214 fBField = f1;
215 fBOn = kTRUE;
216 fNLocal = 5; // helices
217 fRieman = new AliTrackFitterRieman();
218 }
219 else {
220 fBField = 0.0;
221 fBOn = kFALSE;
222 fNLocal = 4;
223 }
224 }
225
483fd7d8 226 if (strstr(st,"SET_PARSIG_TRA")) {
227 sscanf(st,"%s %f",tmp,&f1);
228 fParSigTranslations=f1;
229 }
230
231 if (strstr(st,"SET_PARSIG_ROT")) {
232 sscanf(st,"%s %f",tmp,&f1);
233 fParSigRotations=f1;
234 }
235
236 if (strstr(st,"SET_NSTDDEV")) {
237 sscanf(st,"%s %d",tmp,&idx);
238 fNStdDev=idx;
239 }
240
241 if (strstr(st,"SET_RESCUT_INIT")) {
242 sscanf(st,"%s %f",tmp,&f1);
243 fResCutInitial=f1;
244 }
245
246 if (strstr(st,"SET_RESCUT_OTHER")) {
247 sscanf(st,"%s %f",tmp,&f1);
248 fResCut=f1;
249 }
250
d3603b4d 251 if (strstr(st,"SET_LOCALSIGMAFACTOR")) {
252 sscanf(st,"%s %f %f",tmp,&f1,&f2);
253 if (f1>0 && f2>0) {
254 fSigmaXfactor=f1;
255 fSigmaZfactor=f2;
256 }
257 }
258
483fd7d8 259 if (strstr(st,"SET_STARTFAC")) {
260 sscanf(st,"%s %f",tmp,&f1);
261 fStartFac=f1;
262 }
263
d3603b4d 264 if (strstr(st,"REQUIRE_POINT")) {
265 // syntax: REQUIRE_POINT where ndet updw nreqpts
266 // where = LAYER or DETECTOR
267 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
268 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
269 // nreqpts = minimum number of points of that type
270 sscanf(st,"%s %s %d %d %d",tmp,st2,&itx,&ity,&itz);
271 itx--;
272 if (strstr(st2,"LAYER")) {
273 if (itx<0 || itx>5) return -7;
274 if (ity>0) fNReqLayUp[itx]=itz;
275 else if (ity<0) fNReqLayDown[itx]=itz;
276 else fNReqLay[itx]=itz;
277 fRequirePoints=kTRUE;
278 }
279 else if (strstr(st2,"DETECTOR")) { // DETECTOR
280 if (itx<0 || itx>2) return -7;
281 if (ity>0) fNReqDetUp[itx]=itz;
282 else if (ity<0) fNReqDetDown[itx]=itz;
283 else fNReqDet[itx]=itz;
284 fRequirePoints=kTRUE;
285 }
286 }
287
288
483fd7d8 289 if (strstr(st,"SET_LOCAL_SHIFTS")) { // only enabled mode...
290 fUseLocalShifts = kTRUE;
291 }
292
0856765c 293 if (strstr(st,"MODULE_INDEX")) { // works only for sensitive modules
d3603b4d 294 f1=0; f2=0;
295 sscanf(st,"%s %d %d %d %d %d %d %d %f %f",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips,&f1,&f2);
296 if (idx<0 || idx>2197) return 1; // bad index
483fd7d8 297 voluid=GetModuleVolumeID(idx);
0856765c 298 if (!voluid || voluid>14300) return 1; // bad index
483fd7d8 299 fModuleIndex[nmod]=idx;
300 fModuleVolumeID[nmod]=voluid;
301 fFreeParam[nmod][0]=itx;
302 fFreeParam[nmod][1]=ity;
303 fFreeParam[nmod][2]=itz;
304 fFreeParam[nmod][3]=iph;
305 fFreeParam[nmod][4]=ith;
306 fFreeParam[nmod][5]=ips;
0856765c 307 fMilleModule[nmod] = new AliITSAlignMilleModule(voluid);
d3603b4d 308 if (f1>0) fSensVolSigmaXfactor[idx]=f1;
309 if (f2>0) fSensVolSigmaZfactor[idx]=f2;
483fd7d8 310 nmod++;
311 }
312
313 if (strstr(st,"MODULE_VOLUID")) {
d3603b4d 314 f1=0; f2=0;
315 sscanf(st,"%s %d %d %d %d %d %d %d %f %f",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips,&f1,&f2);
0856765c 316 voluid=UShort_t(idx);
317 if (voluid>14335 && fUseSuperModules) { // custom supermodule
318 int ism=-1;
319 for (int j=0; j<fNSuperModules; j++) {
320 if (voluid==fSuperModule[j]->GetVolumeID()) ism=j;
321 }
322 if (ism<0) return -1; // bad volid
323 fModuleIndex[nmod]=fSuperModule[ism]->GetIndex();
324 fModuleVolumeID[nmod]=voluid;
325 fFreeParam[nmod][0]=itx;
326 fFreeParam[nmod][1]=ity;
327 fFreeParam[nmod][2]=itz;
328 fFreeParam[nmod][3]=iph;
329 fFreeParam[nmod][4]=ith;
330 fFreeParam[nmod][5]=ips;
331 fMilleModule[nmod] = new AliITSAlignMilleModule(*fSuperModule[ism]);
d3603b4d 332 if (f1>0) {
333 for (int kk=0; kk<fMilleModule[nmod]->GetNSensitiveVolumes(); kk++) {
334 idx=AliITSAlignMilleModule::GetIndexFromVolumeID(fMilleModule[nmod]->GetSensitiveVolumeVolumeID()[kk]);
335 if (idx>=0) fSensVolSigmaXfactor[idx]=f1;
336 }
337 }
338 if (f2>0) {
339 for (int kk=0; kk<fMilleModule[nmod]->GetNSensitiveVolumes(); kk++) {
340 idx=AliITSAlignMilleModule::GetIndexFromVolumeID(fMilleModule[nmod]->GetSensitiveVolumeVolumeID()[kk]);
341 if (idx>=0) fSensVolSigmaZfactor[idx]=f2;
342 }
343 } nmod++;
0856765c 344 }
345 else { // sensitive volume
346 idx=GetModuleIndex(voluid);
347 if (idx<0 || idx>2197) return 1; // bad index
348 fModuleIndex[nmod]=idx;
349 fModuleVolumeID[nmod]=voluid;
350 fFreeParam[nmod][0]=itx;
351 fFreeParam[nmod][1]=ity;
352 fFreeParam[nmod][2]=itz;
353 fFreeParam[nmod][3]=iph;
354 fFreeParam[nmod][4]=ith;
355 fFreeParam[nmod][5]=ips;
356 fMilleModule[nmod] = new AliITSAlignMilleModule(voluid);
d3603b4d 357 if (f1>0) fSensVolSigmaXfactor[idx]=f1;
358 if (f2>0) fSensVolSigmaZfactor[idx]=f2;
0856765c 359 nmod++;
360 }
483fd7d8 361 }
0856765c 362 //----------
483fd7d8 363
364 } // end while
365
366 fNModules = nmod;
367 fNGlobal = fNModules*fgNParCh;
368
369 fclose(pfc);
370 return 0;
371}
372
d3603b4d 373void AliITSAlignMille::SetRequiredPoint(Char_t* where, Int_t ndet, Int_t updw, Int_t nreqpts)
374{
375 // set minimum number of points in specific detector or layer
376 // where = LAYER or DETECTOR
377 // ndet = detector number: 1-6 for LAYER and 1-3 for DETECTOR (SPD=1, SDD=2, SSD=3)
378 // updw = 1 for Y>0, -1 for Y<0, 0 if not specified
379 // nreqpts = minimum number of points of that type
380 ndet--;
381 if (strstr(where,"LAYER")) {
382 if (ndet<0 || ndet>5) return;
383 if (updw>0) fNReqLayUp[ndet]=nreqpts;
384 else if (updw<0) fNReqLayDown[ndet]=nreqpts;
385 else fNReqLay[ndet]=nreqpts;
386 fRequirePoints=kTRUE;
387 }
388 else if (strstr(where,"DETECTOR")) {
389 if (ndet<0 || ndet>2) return;
390 if (updw>0) fNReqDetUp[ndet]=nreqpts;
391 else if (updw<0) fNReqDetDown[ndet]=nreqpts;
392 else fNReqDet[ndet]=nreqpts;
393 fRequirePoints=kTRUE;
394 }
395}
396
483fd7d8 397Int_t AliITSAlignMille::GetModuleIndex(const Char_t *symname) {
398 /// index from symname
399 if (!symname) return -1;
400 for (Int_t i=0; i<2198; i++) {
401 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
402 }
403 return -1;
404}
405
406Int_t AliITSAlignMille::GetModuleIndex(UShort_t voluid) {
407 /// index from volume ID
408 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
409 if (lay<1|| lay>6) return -1;
410 Int_t idx=Int_t(voluid)-2048*lay;
411 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
412 for (Int_t ilay=1; ilay<lay; ilay++)
413 idx += AliGeomManager::LayerSize(ilay);
414 return idx;
415}
416
417UShort_t AliITSAlignMille::GetModuleVolumeID(const Char_t *symname) {
418 /// volume ID from symname
0856765c 419 /// works for sensitive volumes only
483fd7d8 420 if (!symname) return 0;
421
422 for (UShort_t voluid=2000; voluid<13300; voluid++) {
423 Int_t modId;
424 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
425 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
426 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
427 }
428 }
429
430 return 0;
431}
432
433UShort_t AliITSAlignMille::GetModuleVolumeID(Int_t index) {
434 /// volume ID from index
0856765c 435 if (index<0) return 0;
436 if (index<2198)
437 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
438 else {
439 for (int i=0; i<fNSuperModules; i++) {
440 if (fSuperModule[i]->GetIndex()==index) return fSuperModule[i]->GetVolumeID();
441 }
442 }
443 return 0;
483fd7d8 444}
445
446void AliITSAlignMille::InitGeometry() {
447 /// initialize geometry
448 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
449 fGeoManager = AliGeomManager::GetGeometry();
450 if (!fGeoManager) {
451 AliInfo("Couldn't initialize geometry");
452 return;
453 }
454 // temporary align object, just use the rotation...
0856765c 455 fTempAlignObj=new AliAlignObjParams;
483fd7d8 456}
457
458void AliITSAlignMille::Init(Int_t nGlobal, /* number of global paramers */
459 Int_t nLocal, /* number of local parameters */
460 Int_t nStdDev /* std dev cut */ )
461{
462 /// Initialization of AliMillepede. Fix parameters, define constraints ...
463 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
464 fIsMilleInit = kTRUE;
465
466 /// Fix non free parameters
467 for (Int_t i=0; i<fNModules; i++) {
468 for (Int_t j=0; j<ITSMILLE_NPARCH; j++) {
469 if (!fFreeParam[i][j]) FixParameter(i*ITSMILLE_NPARCH+j,0.0);
470 else {
d3603b4d 471 // pepopepo: da verificare il settaggio delle sigma, ma forse va bene...
483fd7d8 472 Double_t parsig=0;
473 if (j<3) parsig=fParSigTranslations; // translations (0.0100 cm)
474 else parsig=fParSigRotations; // rotations (1/10 deg)
475 FixParameter(i*ITSMILLE_NPARCH+j,parsig);
476 }
477 }
478 }
479
480
481// // Set iterations
482 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
483}
484
485
486void AliITSAlignMille::AddConstraint(Double_t *par, Double_t value) {
487 /// Constrain equation defined by par to value
488 if (!fIsMilleInit) {
489 AliInfo("Millepede has not been initialized!");
490 return;
491 }
492 fMillepede->SetGlobalConstraint(par, value);
493 AliInfo("Adding constraint");
494}
495
496void AliITSAlignMille::InitGlobalParameters(Double_t *par) {
497 /// Initialize global parameters with par array
498 if (!fIsMilleInit) {
499 AliInfo("Millepede has not been initialized!");
500 return;
501 }
502 fMillepede->SetGlobalParameters(par);
503 AliInfo("Init Global Parameters");
504}
505
506void AliITSAlignMille::FixParameter(Int_t iPar, Double_t value) {
507 /// Parameter iPar is encourage to vary in [-value;value].
508 /// If value == 0, parameter is fixed
509 if (!fIsMilleInit) {
510 AliInfo("Millepede has not been initialized!");
511 return;
512 }
513 fMillepede->SetParSigma(iPar, value);
514 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
515}
516
517void AliITSAlignMille::ResetLocalEquation()
518{
519 /// Reset the derivative vectors
520 for(int i=0; i<fNLocal; i++) {
521 fLocalDerivatives[i] = 0.0;
522 }
523 for(int i=0; i<fNGlobal; i++) {
524 fGlobalDerivatives[i] = 0.0;
525 }
526}
527
0856765c 528Int_t AliITSAlignMille::ApplyToGeometry() {
529 /// apply starting realignment to ideal geometry
530 if (!fGeoManager) return -1;
531 TFile *pref = new TFile(fPreAlignmentFileName.Data());
532 if (!pref->IsOpen()) return -2;
533 TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
534 if (!prea) return -3;
535 Int_t nprea=prea->GetEntriesFast();
536 AliInfo(Form("Array of input misalignments with %d entries",nprea));
537
538 for (int ix=0; ix<nprea; ix++) {
539 AliAlignObjParams *preo=(AliAlignObjParams*) prea->UncheckedAt(ix);
d3603b4d 540 Int_t index=AliITSAlignMilleModule::GetIndexFromVolumeID(preo->GetVolUID());
541 if (index>=0) {
542 fPreAlignQF[index] = (int) preo->GetUniqueID();
543 //printf("index=%d QF=%d\n",index,preo->GetUniqueID());
544 }
0856765c 545 if (!preo->ApplyToGeometry()) return -4;
546 }
547 pref->Close();
548 delete pref;
549
550 fUsePreAlignment = kTRUE;
551 return 0;
552}
553
d3603b4d 554Int_t AliITSAlignMille::GetPreAlignmentQualityFactor(Int_t index) {
555 /// works for sensitive volumes
556 if (!fUsePreAlignment || index<0 || index>2197) return -1;
557 return fPreAlignQF[index];
558}
559
560AliTrackPointArray *AliITSAlignMille::PrepareTrack(AliTrackPointArray *atp) {
561 /// create a new AliTrackPointArray keeping only defined modules
562 /// move points according to a given prealignment, if any
563 /// sort alitrackpoints w.r.t. global Y direction, if selected
564
565 AliTrackPointArray *atps=NULL;
566 Int_t idx[20];
567 Int_t npts=atp->GetNPoints();
568
569 /// checks if AliTrackPoints belong to defined modules
570 Int_t ngoodpts=0;
571 Int_t intidx[20];
572
573 for (int j=0; j<npts; j++) {
574 intidx[j] = IsContained(atp->GetVolumeID()[j]);
575 if (intidx[j]>=0) ngoodpts++;
576 }
577 AliDebug(3,Form("Number of points in defined modules: %d",ngoodpts));
578
579 // reject track if not enough points are left
580 if (ngoodpts<fMinNPtsPerTrack) {
581 AliInfo("Track with not enough points!");
582 return NULL;
583 }
584
585 AliTrackPoint p;
586 // check points in specific places
587 if (fRequirePoints) {
588 Int_t nlayup[6],nlaydown[6],nlay[6];
589 Int_t ndetup[3],ndetdown[3],ndet[3];
590 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
591 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
592
593 for (int i=0; i<npts; i++) {
594 // skip not defined points
595 if (intidx[i]<0) continue;
596 Float_t xx=atp->GetX()[i];
597 Float_t yy=atp->GetY()[i];
598 Float_t r=TMath::Sqrt(xx*xx + yy*yy);
599 int lay=-1;
600 if (r<5) lay=0;
601 else if (r>5 && r<10) lay=1;
602 else if (r>10 && r<18) lay=2;
603 else if (r>18 && r<30) lay=3;
604 else if (r>30 && r<40) lay=4;
605 else if (r>40) lay=5;
606 if (lay<0) continue;
607 int det=lay/2;
608 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
609
610 if (yy>=0.0) { // UP point
611 nlayup[lay]++;
612 nlay[lay]++;
613 ndetup[det]++;
614 ndet[det]++;
615 }
616 else {
617 nlaydown[lay]++;
618 nlay[lay]++;
619 ndetdown[det]++;
620 ndet[det]++;
621 }
622 }
623
624 // checks minimum values
625 Bool_t isok=kTRUE;
626 for (Int_t j=0; j<6; j++) {
627 if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE;
628 if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE;
629 if (nlay[j]<fNReqLay[j]) isok=kFALSE;
630 }
631 for (Int_t j=0; j<3; j++) {
632 if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE;
633 if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE;
634 if (ndet[j]<fNReqDet[j]) isok=kFALSE;
635 }
636 if (!isok) {
637 AliInfo("Track does not meet all location point requirements!");
638 return NULL;
639 }
640 }
641
642 // build a new track with (sorted) (prealigned) good points
643 atps=new AliTrackPointArray(ngoodpts);
644
645 for (int i=0; i<npts; i++) idx[i]=i;
646 // sort track if required
647 if (fUseSortedTracks) TMath::Sort(npts,atp->GetY(),idx); // sort descending...
648
649 Int_t npto=0;
650 for (int i=0; i<npts; i++) {
651 // skip not defined points
652 if (intidx[idx[i]]<0) continue;
653 atp->GetPoint(p,idx[i]);
654
655 // prealign point if required
656 // get IDEAL matrix
657 TGeoHMatrix *svOrigMatrix = fMilleModule[intidx[idx[i]]]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
658 // get back real local coordinates: use OriginalGlobalMatrix because AliTrackPoints were written
659 // with idel geometry
660 Double_t pg[3],pl[3];
661 pg[0]=p.GetX();
662 pg[1]=p.GetY();
663 pg[2]=p.GetZ();
664 AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
665 svOrigMatrix->MasterToLocal(pg,pl);
666
667 AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]));
668
669 // update covariance matrix
670 TGeoHMatrix hcov;
671 Double_t hcovel[9];
672 hcovel[0]=double(p.GetCov()[0]);
673 hcovel[1]=double(p.GetCov()[1]);
674 hcovel[2]=double(p.GetCov()[2]);
675 hcovel[3]=double(p.GetCov()[1]);
676 hcovel[4]=double(p.GetCov()[3]);
677 hcovel[5]=double(p.GetCov()[4]);
678 hcovel[6]=double(p.GetCov()[2]);
679 hcovel[7]=double(p.GetCov()[4]);
680 hcovel[8]=double(p.GetCov()[5]);
681 hcov.SetRotation(hcovel);
682 // now rotate in local system
683 hcov.Multiply(svOrigMatrix);
684 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
685 // now hcov is LOCAL COVARIANCE MATRIX
686
687
688 // pepopepo
689 if (fBug==1) {
690 // correzione bug LAYER 5 SSD temporanea..
691 int ssdidx=AliITSAlignMilleModule::GetIndexFromVolumeID(p.GetVolumeID());
692 if (ssdidx>=500 && ssdidx<1248) {
693 int ladder=(ssdidx-500)%22;
694 if (ladder==18) p.SetVolumeID(AliITSAlignMilleModule::GetVolumeIDFromIndex(ssdidx+1));
695 if (ladder==19) p.SetVolumeID(AliITSAlignMilleModule::GetVolumeIDFromIndex(ssdidx-1));
696 }
697 }
698
699 /// get (evenctually prealigned) matrix of sens. vol.
700 TGeoHMatrix *svMatrix = fMilleModule[intidx[idx[i]]]->GetSensitiveVolumeMatrix(p.GetVolumeID());
701 // modify global coordinates according with pre-aligment
702 svMatrix->LocalToMaster(pl,pg);
703 // now rotate in local system
704 hcov.Multiply(&svMatrix->Inverse());
705 hcov.MultiplyLeft(svMatrix);
706 // hcov is back in GLOBAL RF
707 Float_t pcov[6];
708 pcov[0]=hcov.GetRotationMatrix()[0];
709 pcov[1]=hcov.GetRotationMatrix()[1];
710 pcov[2]=hcov.GetRotationMatrix()[2];
711 pcov[3]=hcov.GetRotationMatrix()[4];
712 pcov[4]=hcov.GetRotationMatrix()[5];
713 pcov[5]=hcov.GetRotationMatrix()[8];
714
715 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
716 AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
717 atps->AddPoint(npto,&p);
718 AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] ));
719
720 npto++;
721 }
722
723 return atps;
724}
725
726
727
728AliTrackPointArray *AliITSAlignMille::SortTrack(AliTrackPointArray *atp) {
729 /// sort alitrackpoints w.r.t. global Y direction
730 AliTrackPointArray *atps=NULL;
731 Int_t idx[20];
732 Int_t npts=atp->GetNPoints();
733 AliTrackPoint p;
734 atps=new AliTrackPointArray(npts);
735
736 TMath::Sort(npts,atp->GetY(),idx);
737
738 for (int i=0; i<npts; i++) {
739 atp->GetPoint(p,idx[i]);
740 atps->AddPoint(i,&p);
741 AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
742 }
743 return atps;
744}
745
746
483fd7d8 747Int_t AliITSAlignMille::InitModuleParams() {
748 /// initialize geometry parameters for a given detector
749 /// for current cluster (fCluster)
750 /// fGlobalInitParam[] is set as:
751 /// [tx,ty,tz,psi,theta,phi]
483fd7d8 752 ///
483fd7d8 753 /// return 0 if success
754
755 if (!fGeoManager) {
756 AliInfo("ITS geometry not initialized!");
757 return -1;
758 }
759
0856765c 760 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
761
483fd7d8 762 // set the internal index (index in module list)
763 UShort_t voluid=fCluster.GetVolumeID();
764 Int_t k=fNModules-1;
d3603b4d 765 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
483fd7d8 766 if (k<0) return -3;
0856765c 767 fCurrentModuleInternalIndex=k; // the internal index of the SUPERMODULE
768
769 fCurrentModuleIndex=fMilleModule[k]->GetIndex(); // index of the SUPERMODULE
483fd7d8 770
771 // set the index
0856765c 772 Int_t index = GetModuleIndex(voluid);
483fd7d8 773 if (index<0) return -2;
0856765c 774 fCurrentSensVolIndex = index; // the index of the SENSITIVE VOLUME
483fd7d8 775
776 fModuleInitParam[0] = 0.0;
777 fModuleInitParam[1] = 0.0;
778 fModuleInitParam[2] = 0.0;
779 fModuleInitParam[3] = 0.0; // psi (X)
780 fModuleInitParam[4] = 0.0; // theta (Y)
781 fModuleInitParam[5] = 0.0; // phi (Z)
0856765c 782
0856765c 783 fCurrentModuleHMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetMatrix();
784
785 for (int ii=0; ii<3; ii++)
786 fCurrentModuleTranslation[ii]=fCurrentModuleHMatrix->GetTranslation()[ii];
787
d3603b4d 788 /// get (evenctually prealigned) matrix of sens. vol.
789 TGeoHMatrix *svMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeMatrix(voluid);
790
483fd7d8 791 fMeasGlo[0] = fCluster.GetX();
792 fMeasGlo[1] = fCluster.GetY();
d3603b4d 793 fMeasGlo[2] = fCluster.GetZ();
794 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
0856765c 795 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
0856765c 796
483fd7d8 797 // set stdev from cluster
798 TGeoHMatrix hcov;
799 Double_t hcovel[9];
800 hcovel[0]=double(fCluster.GetCov()[0]);
801 hcovel[1]=double(fCluster.GetCov()[1]);
d3603b4d 802 hcovel[2]=double(fCluster.GetCov()[2]);
483fd7d8 803 hcovel[3]=double(fCluster.GetCov()[1]);
d3603b4d 804 hcovel[4]=double(fCluster.GetCov()[3]);
483fd7d8 805 hcovel[5]=double(fCluster.GetCov()[4]);
d3603b4d 806 hcovel[6]=double(fCluster.GetCov()[2]);
483fd7d8 807 hcovel[7]=double(fCluster.GetCov()[4]);
808 hcovel[8]=double(fCluster.GetCov()[5]);
809 hcov.SetRotation(hcovel);
810 // now rotate in local system
0856765c 811 hcov.Multiply(svMatrix);
d3603b4d 812 hcov.MultiplyLeft(&svMatrix->Inverse());
483fd7d8 813
d3603b4d 814 // set local sigmas
483fd7d8 815 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
816 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4]));
817 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
818
d3603b4d 819 // set minimum value for SigmaLoc to 10 micron
0856765c 820 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
821 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
822
d3603b4d 823 // multiply local sigmas by global factor
824 fSigmaLoc[0] *= fSigmaXfactor;
825 fSigmaLoc[2] *= fSigmaZfactor;
826
827 // multiply local sigmas by individual factor
828 fSigmaLoc[0] *= fSensVolSigmaXfactor[index];
829 fSigmaLoc[2] *= fSensVolSigmaZfactor[index];
830
831 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
483fd7d8 832
833 return 0;
834}
835
836void AliITSAlignMille::SetCurrentModule(Int_t index) {
0856765c 837 /// set as current the SuperModule that contains the 'index' sens.vol.
838 if (index<0 || index>2197) {
839 AliInfo("index does not correspond to a sensitive volume!");
840 return;
841 }
483fd7d8 842 UShort_t voluid=GetModuleVolumeID(index);
0856765c 843 Int_t k=IsContained(voluid);
844 if (k>=0){
0856765c 845 fCluster.SetVolumeID(voluid);
0856765c 846 fCluster.SetXYZ(0,0,0);
847 InitModuleParams();
848 }
849 else
d3603b4d 850 AliInfo(Form("module %d not defined\n",index));
0856765c 851}
852
853void AliITSAlignMille::SetCurrentSensitiveModule(Int_t index) {
854 /// set as current the SuperModule that contains the 'index' sens.vol.
855 if (index<0 || index>2197) {
856 AliInfo("index does not correspond to a sensitive volume!");
857 return;
858 }
859 UShort_t voluid=AliITSAlignMilleModule::GetVolumeIDFromIndex(index);
860 Int_t k=IsDefined(voluid);
0856765c 861 if (k>=0){
483fd7d8 862 fCluster.SetVolumeID(voluid);
863 fCluster.SetXYZ(0,0,0);
864 InitModuleParams();
865 }
0856765c 866 else
d3603b4d 867 AliInfo(Form("module %d not defined\n",index));
483fd7d8 868}
869
0856765c 870void AliITSAlignMille::Print(Option_t*) const
871{
d3603b4d 872 /// print infos
483fd7d8 873 printf("*** AliMillepede for ITS ***\n");
0856765c 874 printf(" number of defined super modules: %d\n",fNModules);
875
483fd7d8 876 if (fGeoManager)
877 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
878 else
879 printf(" geometry not loaded\n");
0856765c 880
881 if (fUseSuperModules)
882 printf(" using custom supermodules ( %d defined )\n",fNSuperModules);
883 else
884 printf(" custom supermodules not used\n");
885
886 if (fUsePreAlignment)
887 printf(" using prealignment from %s \n",fPreAlignmentFileName.Data());
888 else
889 printf(" prealignment not used\n");
890
d3603b4d 891 if (fBOn)
892 printf(" B Field set to %f T - using Riemann's helices\n",fBField);
893 else
894 printf(" B Field OFF - using straight lines \n");
895
483fd7d8 896 if (fUseLocalShifts)
897 printf(" Alignment shifts will be computed in LOCAL RS\n");
898 else
899 printf(" Alignment shifts will be computed in GLOBAL RS\n");
d3603b4d 900
901 if (fRequirePoints) printf(" Required points in tracks:\n");
902 for (Int_t i=0; i<6; i++) {
903 if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
904 if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
905 if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]);
906 }
907 for (Int_t i=0; i<3; i++) {
908 if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
909 if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
910 if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]);
911 }
0856765c 912
d3603b4d 913 printf("\n Millepede configuration parameters:\n");
914 printf(" init parsig for translations : %.4f\n",fParSigTranslations);
915 printf(" init parsig for rotations : %.4f\n",fParSigRotations);
916 printf(" init value for chi2 cut : %.4f\n",fStartFac);
917 printf(" first iteration cut value : %.4f\n",fResCutInitial);
918 printf(" other iterations cut value : %.4f\n",fResCut);
919 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
920 printf(" mult. fact. for local sigmas : %.4f %.4f\n",fSigmaXfactor, fSigmaZfactor);
483fd7d8 921
922 printf("List of defined modules:\n");
0856765c 923 printf(" intidx\tindex\tvoluid\tname\n");
483fd7d8 924 for (int i=0; i<fNModules; i++)
0856765c 925 printf(" %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],fMilleModule[i]->GetName());
926
483fd7d8 927}
928
0856765c 929AliITSAlignMilleModule *AliITSAlignMille::GetMilleModule(UShort_t voluid)
930{
931 // return pointer to a define supermodule
932 // return NULL if error
933 Int_t i=IsDefined(voluid);
934 if (i<0) return NULL;
935 return fMilleModule[i];
936}
937
938AliITSAlignMilleModule *AliITSAlignMille::GetCurrentModule()
939{
940 if (fNModules) return fMilleModule[fCurrentModuleInternalIndex];
941 return NULL;
942}
943
944void AliITSAlignMille::PrintCurrentModuleInfo()
945{
483fd7d8 946 ///
0856765c 947 Int_t k=fCurrentModuleInternalIndex;
948 if (k<0 || k>=fNModules) return;
949 fMilleModule[k]->Print("");
483fd7d8 950}
951
d3603b4d 952Bool_t AliITSAlignMille::InitRiemanFit() {
953 // Initialize Riemann Fitter for current track
954 // return kFALSE if error
955
956 if (!fBOn) return kFALSE;
957
958 Int_t npts=0;
959 AliTrackPoint ap;
960 npts = fTrack->GetNPoints();
961 AliDebug(3,Form("Fitting track with %d points",npts));
962
963 fRieman->Reset();
964 fRieman->SetTrackPointArray(fTrack);
965
966 TArrayI ai(npts);
967 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
968
969 // fit track with 5 params in his own tracking-rotated reference system
970 // xc = -p[1]/p[0];
971 // yc = 1/p[0];
972 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
973 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
974 return kFALSE;
975 }
976
977 for (int i=0; i<5; i++)
978 fLocalInitParam[i] = fRieman->GetParam()[i];
979
980 return kTRUE;
981}
982
483fd7d8 983
984void AliITSAlignMille::InitTrackParams(int meth) {
985 /// initialize local parameters with different methods
986 /// for current track (fTrack)
987
0856765c 988 Int_t npts=0;
989 TF1 *f1=NULL;
990 TGraph *g=NULL;
991 Float_t sigmax[20],sigmay[20],sigmaz[20];
992 AliTrackPoint ap;
993 TGraphErrors *ge=NULL;
994
483fd7d8 995 switch (meth) {
996 case 1: // simple linear interpolation
997 // get local starting parameters (to be substituted by ESD track parms)
998 // local parms (fLocalInitParam[]) are:
999 // [0] = global x coord. of straight line intersection at y=0 plane
1000 // [1] = global z coord. of straight line intersection at y=0 plane
1001 // [2] = px/py
1002 // [3] = pz/py
1003
1004 // test #1: linear fit in x(y) and z(y)
0856765c 1005 npts = fTrack->GetNPoints();
d3603b4d 1006 AliDebug(3,Form("*** initializing track with %d points ***",npts));
483fd7d8 1007
0856765c 1008 f1=new TF1("f1","[0]+x*[1]",-50,50);
483fd7d8 1009
0856765c 1010 g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
483fd7d8 1011 g->Fit(f1,"RNQ");
1012 fLocalInitParam[0] = f1->GetParameter(0);
1013 fLocalInitParam[2] = f1->GetParameter(1);
1014 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
1015 delete g; g=NULL;
1016
1017 g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ());
1018 g->Fit(f1,"RNQ");
1019 fLocalInitParam[1] = f1->GetParameter(0);
1020 fLocalInitParam[3] = f1->GetParameter(1);
1021 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
1022 delete g;
1023 delete f1;
1024
1025 break;
0856765c 1026
1027 case 2: // simple linear interpolation weighted using sigmas
1028 // get local starting parameters (to be substituted by ESD track parms)
1029 // local parms (fLocalInitParam[]) are:
1030 // [0] = global x coord. of straight line intersection at y=0 plane
1031 // [1] = global z coord. of straight line intersection at y=0 plane
1032 // [2] = px/py
1033 // [3] = pz/py
1034
1035 // test #1: linear fit in x(y) and z(y)
1036 npts = fTrack->GetNPoints();
d3603b4d 1037 AliDebug(3,Form("*** initializing track with %d points ***",npts));
1038
0856765c 1039 for (Int_t isig=0; isig<npts; isig++) {
1040 fTrack->GetPoint(ap,isig);
1041 sigmax[isig]=ap.GetCov()[0];
1042 if (sigmax[isig]<1.0e-07) sigmax[isig]=1.0e-07; // minimum sigma=3 mu
1043 sigmax[isig]=TMath::Sqrt(sigmax[isig]);
1044
d3603b4d 1045 sigmay[isig]=ap.GetCov()[3];
0856765c 1046 if (sigmay[isig]<1.0e-07) sigmay[isig]=1.0e-07; // minimum sigma=3 mu
1047 sigmay[isig]=TMath::Sqrt(sigmay[isig]);
1048
1049 sigmaz[isig]=ap.GetCov()[5];
1050 if (sigmaz[isig]<1.0e-07) sigmaz[isig]=1.0e-07; // minimum sigma=3 mu
1051 sigmaz[isig]=TMath::Sqrt(sigmaz[isig]);
d3603b4d 1052
1053 if (fTempExcludedModule>=0 && fTempExcludedModule==AliITSAlignMilleModule::GetIndexFromVolumeID(ap.GetVolumeID())) {
1054 sigmax[isig] *= 1000.;
1055 sigmay[isig] *= 1000.;
1056 sigmaz[isig] *= 1000.;
1057 fTempExcludedModule=-1;
1058 }
0856765c 1059 }
1060
1061 f1=new TF1("f1","[0]+x*[1]",-50,50);
1062
1063 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetX(),sigmay,sigmax);
1064 ge->Fit(f1,"RNQ");
1065 fLocalInitParam[0] = f1->GetParameter(0);
1066 fLocalInitParam[2] = f1->GetParameter(1);
1067 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
1068 delete ge; ge=NULL;
1069
1070 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetZ(),sigmay,sigmaz);
1071 ge->Fit(f1,"RNQ");
1072 fLocalInitParam[1] = f1->GetParameter(0);
1073 fLocalInitParam[3] = f1->GetParameter(1);
1074 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
1075 delete ge;
1076 delete f1;
1077
1078 break;
1079
483fd7d8 1080 }
0856765c 1081}
483fd7d8 1082
0856765c 1083Int_t AliITSAlignMille::IsDefined(UShort_t voluid) const
1084{
1085 // checks if supermodule 'voluid' is defined and return the internal index
1086 // return -1 if error
1087 Int_t k=fNModules-1;
1088 while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
1089 if (k<0) return -1;
1090 return k;
1091}
1092
1093Int_t AliITSAlignMille::IsContained(UShort_t voluid) const
1094{
1095 // checks if the sensitive module 'voluid' is contained inside a supermodule and return the internal index of the last identified supermodule
1096 // return -1 if error
1097 if (AliITSAlignMilleModule::GetIndexFromVolumeID(voluid)<0) return -1;
1098 Int_t k=fNModules-1;
1099 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
1100 if (k<0) return -1;
1101 return k;
483fd7d8 1102}
0856765c 1103
483fd7d8 1104Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
1105{
0856765c 1106 /// check if a sensitive volume is contained inside one of the defined supermodules
483fd7d8 1107 Int_t k=fNModules-1;
0856765c 1108 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
483fd7d8 1109 if (k>=0) return kTRUE;
1110 return kFALSE;
1111}
1112
1113Int_t AliITSAlignMille::CheckCurrentTrack() {
1114 /// checks if AliTrackPoints belongs to defined modules
1115 /// return number of good poins
1116 /// return 0 if not enough points
1117
1118 Int_t npts = fTrack->GetNPoints();
1119 Int_t ngoodpts=0;
1120 // debug points
1121 for (int j=0; j<npts; j++) {
1122 UShort_t voluid = fTrack->GetVolumeID()[j];
1123 if (CheckVolumeID(voluid)) {
1124 ngoodpts++;
1125 }
1126 }
d3603b4d 1127
483fd7d8 1128 if (ngoodpts<fMinNPtsPerTrack) return 0;
1129
1130 return ngoodpts;
1131}
1132
1133Int_t AliITSAlignMille::ProcessTrack(AliTrackPointArray *track) {
1134 /// Process track; Loop over hits and set local equations
1135 /// here 'track' is a AliTrackPointArray
1136 /// return 0 if success;
1137
1138 if (!fIsMilleInit) {
1139 AliInfo("Millepede has not been initialized!");
1140 return -1;
1141 }
1142
1143 Int_t npts = track->GetNPoints();
d3603b4d 1144 AliDebug(2,Form("*** Input track with %d points ***",npts));
483fd7d8 1145
d3603b4d 1146 // preprocessing of the input track: keep only points in defined volumes,
1147 // move points if prealignment is set, sort by Yglo if required
1148
1149 fTrack=PrepareTrack(track);
1150 if (!fTrack) return -1;
483fd7d8 1151
d3603b4d 1152 npts = fTrack->GetNPoints();
1153 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
1154
1155 if (!fBOn) { // straight lines
1156 // set local starting parameters (to be substituted by ESD track parms)
1157 // local parms (fLocalInitParam[]) are:
1158 // [0] = global x coord. of straight line intersection at y=0 plane
1159 // [1] = global z coord. of straight line intersection at y=0 plane
1160 // [2] = px/py
1161 // [3] = pz/py
1162 InitTrackParams(fInitTrackParamsMeth);
1163 }
1164 else {
1165 // local parms (fLocalInitParam[]) are the Riemann Fitter params
1166 if (!InitRiemanFit()) {
1167 AliInfo("Riemann fit failed! skipping this track...");
1168 delete fTrack;
1169 fTrack=NULL;
1170 return -5;
1171 }
1172 }
1173
1174 Int_t nloceq=0;
1175 MilleData *md = new MilleData[npts];
1176
483fd7d8 1177 for (Int_t ipt=0; ipt<npts; ipt++) {
1178 fTrack->GetPoint(fCluster,ipt);
d3603b4d 1179 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
483fd7d8 1180
1181 // set geometry parameters for the the current module
483fd7d8 1182 if (InitModuleParams()) continue;
483fd7d8 1183 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
d3603b4d 1184 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
1185
1186 if (!AddLocalEquation(md[nloceq])) {
1187 nloceq++;
1188 fProcessedPoints[fCurrentModuleInternalIndex]++;
1189 }
1190 else {
1191 fTotBadLocEqPoints++;
1192 }
483fd7d8 1193
483fd7d8 1194 } // end loop over points
1195
d3603b4d 1196 delete fTrack;
1197 fTrack=NULL;
1198
1199 // not enough good points!
1200 if (nloceq<fMinNPtsPerTrack) {
1201 delete [] md;
1202 return -1;
1203 }
1204
1205 // finally send local equations to millepede
1206 SetLocalEquations(md,nloceq);
1207
1208 delete [] md;
1209
483fd7d8 1210 return 0;
1211}
1212
1213Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) {
d3603b4d 1214 /// calculate intersection point of track with current module in local coordinates
1215 /// according with a given set of parameters (local(4/5) and global(6))
483fd7d8 1216 /// and fill fPintLoc/Glo
d3603b4d 1217 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
483fd7d8 1218 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
1219 /// return 0 if success
1220
d3603b4d 1221 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]));
1222 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]));
483fd7d8 1223
483fd7d8 1224
d3603b4d 1225 // prepare the TGeoHMatrix
0856765c 1226 TGeoHMatrix *fTempHMat = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar);
1227 if (!fTempHMat) return -1;
d3603b4d 1228
1229 Double_t v0g[3]; // vector with straight line direction in global coord.
1230 Double_t p0g[3]; // point of the straight line (glo)
1231
1232 if (fBOn) { // B FIELD!
1233 AliTrackPoint prf;
1234 for (int ip=0; ip<5; ip++)
1235 fRieman->SetParam(ip,lpar[ip]);
1236
1237 if (!fRieman->GetPCA(fCluster,prf)) {
1238 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
1239 return -3;
1240 }
1241 // now determine straight line passing tangent to fit curve at prf
1242 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
1243 // mo' P1=(X,Y,Z)_glo_prf
1244 // => (x,y,Z)_trk_prf ruotando di alpha...
1245 Double_t alpha=fRieman->GetAlpha();
1246 Double_t x1g = prf.GetX();
1247 Double_t y1g = prf.GetY();
1248 Double_t z1g = prf.GetZ();
1249 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
1250 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
1251 Double_t z1t = z1g;
1252
1253 Double_t x2t = x1t+1.0;
1254 Double_t y2t = y1t+fRieman->GetDYat(x1t);
1255 Double_t z2t = z1t+fRieman->GetDZat(x1t);
1256 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
1257 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
1258 Double_t z2g = z2t;
1259
1260 AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi()));
1261 AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
1262 AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
1263
1264 if (TMath::Abs(y2g-y1g)<1e-15) {
1265 AliInfo("DeltaY=0! Cannot proceed...");
1266 return -1;
1267 }
1268 // ugx,1,ugz
1269 v0g[0] = (x2g-x1g)/(y2g-y1g);
1270 v0g[1]=1.0;
1271 v0g[2] = (z2g-z1g)/(y2g-y1g);
1272
1273 // point: just keep prf
1274 p0g[0]=x1g;
1275 p0g[1]=y1g;
1276 p0g[2]=z1g;
1277 }
1278 else { // staight line
1279 // vector of initial straight line direction in glob. coord
1280 v0g[0]=lpar[2];
1281 v0g[1]=1.0;
1282 v0g[2]=lpar[3];
1283
1284 // intercept in yg=0 plane in glob coord
1285 p0g[0]=lpar[0];
1286 p0g[1]=0.0;
1287 p0g[2]=lpar[1];
1288 }
1289 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]));
1290
483fd7d8 1291 // same in local coord.
1292 Double_t p0l[3],v0l[3];
1293 fTempHMat->MasterToLocalVect(v0g,v0l);
1294 fTempHMat->MasterToLocal(p0g,p0l);
d3603b4d 1295
483fd7d8 1296 if (TMath::Abs(v0l[1])<1e-15) {
1297 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
1298 return -1;
1299 }
1300
1301 // local intersection point
1302 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
1303 fPintLoc[1] = 0;
1304 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
1305
1306 // global intersection point
1307 fTempHMat->LocalToMaster(fPintLoc,fPintGlo);
1308 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]));
d3603b4d 1309
483fd7d8 1310 return 0;
1311}
1312
1313Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) {
1314 /// calculate numerically (ROOT's style) the derivatives for
1315 /// local X intersection and local Z intersection
d3603b4d 1316 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
483fd7d8 1317 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
1318 /// return 0 if success
1319
1320 // copy initial parameters
1321 Double_t lpar[ITSMILLE_NLOCAL];
1322 Double_t gpar[ITSMILLE_NPARCH];
1323 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) lpar[i]=fLocalInitParam[i];
1324 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) gpar[i]=fModuleInitParam[i];
1325
483fd7d8 1326 // trial with fixed dpar...
1327 Double_t dpar=0.0;
d3603b4d 1328
1329 if (islpar) { // track parameters
483fd7d8 1330 //dpar=fLocalInitParam[paridx]*0.001;
1331 // set minimum dpar
d3603b4d 1332 if (!fBOn) {
1333 if (paridx<2) dpar=1.0e-4; // translations
1334 else dpar=1.0e-6; // direction
1335 }
1336 else { // B Field
1337 // pepo: proviamo con 1/1000, poi evenctually 1/100...
1338 Double_t dfrac=0.01;
1339 switch(paridx) {
1340 case 0:
1341 // RMS cosmics: 1e-4
1342 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1343 break;
1344 case 1:
1345 // RMS cosmics: 0.2
1346 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1347 break;
1348 case 2:
1349 // RMS cosmics: 9
1350 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1351 break;
1352 case 3:
1353 // RMS cosmics: 7
1354 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1355 break;
1356 case 4:
1357 // RMS cosmics: 0.3
1358 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1359 break;
1360 }
1361 }
483fd7d8 1362 }
d3603b4d 1363 else { // alignment global parameters
483fd7d8 1364 //dpar=fModuleInitParam[paridx]*0.001;
1365 if (paridx<3) dpar=1.0e-4; // translations
1366 else dpar=1.0e-2; // angles
1367 }
d3603b4d 1368
1369 AliDebug(3,Form("+++ using dpar=%g",dpar));
483fd7d8 1370
1371 // calculate derivative ROOT's like:
1372 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1373 Double_t pintl1[3]; // f(x-h)
1374 Double_t pintl2[3]; // f(x-h/2)
1375 Double_t pintl3[3]; // f(x+h/2)
1376 Double_t pintl4[3]; // f(x+h)
1377
1378 // first values
1379 if (islpar) lpar[paridx] -= dpar;
1380 else gpar[paridx] -= dpar;
1381 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1382 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
1383
1384 // second values
1385 if (islpar) lpar[paridx] += dpar/2;
1386 else gpar[paridx] += dpar/2;
1387 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1388 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
1389
1390 // third values
1391 if (islpar) lpar[paridx] += dpar;
1392 else gpar[paridx] += dpar;
1393 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1394 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
1395
1396 // fourth values
1397 if (islpar) lpar[paridx] += dpar/2;
1398 else gpar[paridx] += dpar/2;
1399 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1400 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
1401
1402 Double_t h2 = 1./(2.*dpar);
1403 Double_t d0 = pintl4[0]-pintl1[0];
1404 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
1405 fDerivativeXLoc = h2*(4*d2 - d0)/3.;
1406 if (TMath::Abs(fDerivativeXLoc) < 1.0e-9) fDerivativeXLoc=0.0;
1407
1408 d0 = pintl4[2]-pintl1[2];
1409 d2 = 2.*(pintl3[2]-pintl2[2]);
1410 fDerivativeZLoc = h2*(4*d2 - d0)/3.;
1411 if (TMath::Abs(fDerivativeZLoc) < 1.0e-9) fDerivativeZLoc=0.0;
1412
1413 AliDebug(3,Form("\n+++ derivatives +++ \n"));
1414 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",fDerivativeXLoc));
1415 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",fDerivativeZLoc));
1416
1417 return 0;
1418}
1419
d3603b4d 1420
1421Int_t AliITSAlignMille::AddLocalEquation(MilleData &m) {
1422 /// Define local equation for current cluster in X and Z coor.
1423 /// and store them to memory
483fd7d8 1424 /// return 0 if success
1425
1426 // store first interaction point
d3603b4d 1427 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -4;
483fd7d8 1428 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
0856765c 1429 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
483fd7d8 1430
1431 // calculate local derivatives numerically
1432 Double_t dXdL[ITSMILLE_NLOCAL],dZdL[ITSMILLE_NLOCAL];
d3603b4d 1433 for (Int_t i=0; i<fNLocal; i++) {
483fd7d8 1434 if (CalcDerivatives(i,kTRUE)) return -1;
1435 dXdL[i]=fDerivativeXLoc;
1436 dZdL[i]=fDerivativeZLoc;
1437 }
1438
1439 Double_t dXdG[ITSMILLE_NPARCH],dZdG[ITSMILLE_NPARCH];
1440 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
1441 if (CalcDerivatives(i,kFALSE)) return -1;
1442 dXdG[i]=fDerivativeXLoc;
1443 dZdG[i]=fDerivativeZLoc;
1444 }
1445
1446 AliDebug(2,Form("\n***************\n"));
d3603b4d 1447 for (Int_t i=0; i<fNLocal; i++)
483fd7d8 1448 AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
1449 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1450 AliDebug(2,Form("Global parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdG[i],dZdG[i]));
1451 AliDebug(2,Form("\n***************\n"));
1452
d3603b4d 1453 // check if at least 1 local and 1 global derivs. are not null
1454 Double_t nonzero=0.0;
1455 for (Int_t i=0; i<fNLocal; i++) nonzero += TMath::Abs(dXdL[i]);
1456 if (nonzero==0.0) {
1457 AliInfo("Aborting local equations for this point beacuse of zero local X derivatives!");
1458 return -2;
1459 }
1460 nonzero=0.0;
1461 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) nonzero += TMath::Abs(dXdG[i]);
1462 if (nonzero==0.0) {
1463 AliInfo("Aborting local equations for this point beacuse of zero global X derivatives!");
1464 return -2;
1465 }
1466 nonzero=0.0;
1467 for (Int_t i=0; i<fNLocal; i++) nonzero += TMath::Abs(dZdL[i]);
1468 if (nonzero==0.0) {
1469 AliInfo("Aborting local equations for this point beacuse of zero local Z derivatives!");
1470 return -2;
1471 }
1472 nonzero=0.0;
1473 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) nonzero += TMath::Abs(dZdG[i]);
1474 if (nonzero==0.0) {
1475 AliInfo("Aborting local equations for this point beacuse of zero global Z derivatives!");
1476 return -2;
1477 }
1478
1479 // ok, can copy to m
1480
1481 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
483fd7d8 1482 // set equation for Xloc coordinate
d3603b4d 1483 for (Int_t i=0; i<fNLocal; i++) {
1484 m.idxlocX[i]=i;
1485 m.derlocX[i]=dXdL[i];
1486 }
1487 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
1488 m.idxgloX[i]=fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i;
1489 m.dergloX[i]=dXdG[i];
1490 }
1491 m.measX = fMeasLoc[0]-fPintLoc0[0];
1492 m.sigmaX = fSigmaLoc[0];
483fd7d8 1493
d3603b4d 1494 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
483fd7d8 1495 // set equation for Zloc coordinate
d3603b4d 1496 for (Int_t i=0; i<fNLocal; i++) {
1497 m.idxlocZ[i]=i;
1498 m.derlocZ[i]=dZdL[i];
1499 }
1500 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
1501 m.idxgloZ[i]=fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i;
1502 m.dergloZ[i]=dZdG[i];
1503 }
1504 m.measZ = fMeasLoc[2]-fPintLoc0[2];
1505 m.sigmaZ = fSigmaLoc[2];
1506
483fd7d8 1507 return 0;
1508}
1509
1510
d3603b4d 1511void AliITSAlignMille::SetLocalEquations(MilleData *m, Int_t neq) {
1512 /// Set local equations with data stored in m
1513 /// return 0 if success
1514
1515 for (Int_t j=0; j<neq; j++) {
1516
1517 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",m[j].measX, m[j].sigmaX));
1518 // set equation for Xloc coordinate
1519 for (Int_t i=0; i<fNLocal; i++)
1520 SetLocalDerivative( m[j].idxlocX[i], m[j].derlocX[i] );
1521
1522 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1523 SetGlobalDerivative( m[j].idxgloX[i], m[j].dergloX[i] );
1524
1525 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m[j].measX, m[j].sigmaX);
1526
1527
1528 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",m[j].measZ, m[j].sigmaZ));
1529 // set equation for Zloc coordinate
1530 for (Int_t i=0; i<fNLocal; i++)
1531 SetLocalDerivative( m[j].idxlocZ[i], m[j].derlocZ[i] );
1532
1533 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
1534 SetGlobalDerivative( m[j].idxgloZ[i], m[j].dergloZ[i] );
1535
1536 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m[j].measZ, m[j].sigmaZ);
1537 }
1538}
1539
1540
483fd7d8 1541void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
1542 /// Call local fit for this track
1543 if (!fIsMilleInit) {
1544 AliInfo("Millepede has not been initialized!");
1545 return;
1546 }
1547 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
1548 AliDebug(2,Form("iRes = %d",iRes));
1549 if (iRes && !lSingleFit) {
1550 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
1551 }
1552}
1553
1554void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
1555 /// Call global fit; Global parameters are stored in parameters
1556 if (!fIsMilleInit) {
1557 AliInfo("Millepede has not been initialized!");
1558 return;
1559 }
1560 fMillepede->GlobalFit(parameters,errors,pulls);
1561 AliInfo("Done fitting global parameters!");
1562}
1563
1564Double_t AliITSAlignMille::GetParError(Int_t iPar) {
1565 /// Get error of parameter iPar
1566 if (!fIsMilleInit) {
1567 AliInfo("Millepede has not been initialized!");
1568 return 0;
1569 }
1570 Double_t lErr = fMillepede->GetParError(iPar);
1571 return lErr;
1572}
1573
1574void AliITSAlignMille::PrintGlobalParameters() {
1575 /// Print global parameters
1576 if (!fIsMilleInit) {
1577 AliInfo("Millepede has not been initialized!");
1578 return;
1579 }
1580 fMillepede->PrintGlobalParameters();
1581}
1582
1583// //_________________________________________________________________________
0856765c 1584Int_t AliITSAlignMille::LoadSuperModuleFile(const Char_t *sfile)
1585{
1586 // load definitions of supermodules from a root file
1587 // return 0 if success
1588
1589 TFile *smf=TFile::Open(sfile);
1590 if (!smf->IsOpen()) {
1591 AliInfo(Form("Cannot open supermodule file %s",sfile));
1592 return -1;
1593 }
1594
1595 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
1596 if (!sma) {
1597 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
1598 return -2;
1599 }
1600 Int_t nsma=sma->GetEntriesFast();
1601 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
1602
1603 Char_t st[250];
1604 char symname[150];
1605 UShort_t volid;
1606 TGeoHMatrix m;
1607
1608 for (Int_t i=0; i<nsma; i++) {
1609 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
1610 volid=a->GetVolUID();
1611 strcpy(st,a->GetSymName());
1612 a->GetMatrix(m);
1613
1614 sscanf(st,"%s",symname);
1615 // decode module list
1616 char *stp=strstr(st,"ModuleList:");
1617 if (!stp) return -3;
1618 stp += 11;
1619 int idx[2200];
1620 char spp[200]; int jp=0;
1621 char cl[20];
1622 strcpy(st,stp);
1623 int l=strlen(st);
1624 int j=0;
1625 int n=0;
1626
1627 while (j<=l) {
1628 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
1629 spp[jp]=0;
1630 jp=0;
1631 if (strlen(spp)) {
1632 int k=strcspn(spp,"-");
1633 if (k<int(strlen(spp))) { // c'e' il -
1634 strcpy(cl,&(spp[k+1]));
1635 spp[k]=0;
1636 int ifrom=atoi(spp); int ito=atoi(cl);
1637 for (int b=ifrom; b<=ito; b++) {
1638 idx[n]=b;
1639 n++;
1640 }
1641 }
1642 else { // numerillo singolo
1643 idx[n]=atoi(spp);
1644 n++;
1645 }
1646 }
1647 }
1648 else {
1649 spp[jp]=st[j];
1650 jp++;
1651 }
1652 j++;
1653 }
1654 UShort_t volidsv[2198];
1655 for (j=0;j<n;j++) {
1656 volidsv[j]=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx[j]);
1657 if (!volidsv[j]) {
1658 AliInfo(Form("Index %d not valid (range 0->2197)",idx[j]));
1659 return -5;
1660 }
1661 }
1662 Int_t smindex=int(2198+volid-14336); // virtual index
1663 fSuperModule[fNSuperModules]=new AliITSAlignMilleModule(smindex,volid,symname,&m,n,volidsv);
1664
1665 //-------------
1666 fNSuperModules++;
1667 }
1668
1669 smf->Close();
1670
1671 fUseSuperModules=1;
1672 return 0;
1673}
483fd7d8 1674