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