]> 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
Fixes for bug #61900: Restore capaibility to chooses ITS PID method (Marco)
[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")) {
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++) {
75d480f6 468 for (Int_t j=0; j<ITSMILLENPARCH; j++) {
469 if (!fFreeParam[i][j]) FixParameter(i*ITSMILLENPARCH+j,0.0);
483fd7d8 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)
75d480f6 475 FixParameter(i*ITSMILLENPARCH+j,parsig);
483fd7d8 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
75d480f6 528// newpep
0856765c 529Int_t AliITSAlignMille::ApplyToGeometry() {
530 /// apply starting realignment to ideal geometry
75d480f6 531 if(!AliGeomManager::GetGeometry()) return -1;
532
0856765c 533 TFile *pref = new TFile(fPreAlignmentFileName.Data());
534 if (!pref->IsOpen()) return -2;
535 TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
536 if (!prea) return -3;
537 Int_t nprea=prea->GetEntriesFast();
538 AliInfo(Form("Array of input misalignments with %d entries",nprea));
539
75d480f6 540 AliGeomManager::ApplyAlignObjsToGeom(*prea); // apply all levels of objs
541
542 // set prealignment factor if defined...
0856765c 543 for (int ix=0; ix<nprea; ix++) {
544 AliAlignObjParams *preo=(AliAlignObjParams*) prea->UncheckedAt(ix);
d3603b4d 545 Int_t index=AliITSAlignMilleModule::GetIndexFromVolumeID(preo->GetVolUID());
546 if (index>=0) {
547 fPreAlignQF[index] = (int) preo->GetUniqueID();
548 //printf("index=%d QF=%d\n",index,preo->GetUniqueID());
549 }
75d480f6 550 //if (!preo->ApplyToGeometry()) return -4;
0856765c 551 }
552 pref->Close();
553 delete pref;
554
555 fUsePreAlignment = kTRUE;
556 return 0;
557}
75d480f6 558// endnewpep
0856765c 559
75d480f6 560Int_t AliITSAlignMille::GetPreAlignmentQualityFactor(Int_t index) const {
d3603b4d 561 /// works for sensitive volumes
562 if (!fUsePreAlignment || index<0 || index>2197) return -1;
563 return fPreAlignQF[index];
564}
565
566AliTrackPointArray *AliITSAlignMille::PrepareTrack(AliTrackPointArray *atp) {
567 /// create a new AliTrackPointArray keeping only defined modules
568 /// move points according to a given prealignment, if any
569 /// sort alitrackpoints w.r.t. global Y direction, if selected
570
571 AliTrackPointArray *atps=NULL;
572 Int_t idx[20];
573 Int_t npts=atp->GetNPoints();
574
575 /// checks if AliTrackPoints belong to defined modules
576 Int_t ngoodpts=0;
577 Int_t intidx[20];
578
579 for (int j=0; j<npts; j++) {
580 intidx[j] = IsContained(atp->GetVolumeID()[j]);
581 if (intidx[j]>=0) ngoodpts++;
582 }
583 AliDebug(3,Form("Number of points in defined modules: %d",ngoodpts));
584
585 // reject track if not enough points are left
586 if (ngoodpts<fMinNPtsPerTrack) {
587 AliInfo("Track with not enough points!");
588 return NULL;
589 }
590
591 AliTrackPoint p;
592 // check points in specific places
593 if (fRequirePoints) {
594 Int_t nlayup[6],nlaydown[6],nlay[6];
595 Int_t ndetup[3],ndetdown[3],ndet[3];
596 for (Int_t j=0; j<6; j++) {nlayup[j]=0; nlaydown[j]=0; nlay[j]=0;}
597 for (Int_t j=0; j<3; j++) {ndetup[j]=0; ndetdown[j]=0; ndet[j]=0;}
598
599 for (int i=0; i<npts; i++) {
600 // skip not defined points
601 if (intidx[i]<0) continue;
602 Float_t xx=atp->GetX()[i];
603 Float_t yy=atp->GetY()[i];
604 Float_t r=TMath::Sqrt(xx*xx + yy*yy);
605 int lay=-1;
606 if (r<5) lay=0;
607 else if (r>5 && r<10) lay=1;
608 else if (r>10 && r<18) lay=2;
609 else if (r>18 && r<30) lay=3;
610 else if (r>30 && r<40) lay=4;
611 else if (r>40) lay=5;
612 if (lay<0) continue;
613 int det=lay/2;
614 //printf("Point %d - x=%f y=%f R=%f lay=%d det=%d\n",i,xx,yy,r,lay,det);
615
616 if (yy>=0.0) { // UP point
617 nlayup[lay]++;
618 nlay[lay]++;
619 ndetup[det]++;
620 ndet[det]++;
621 }
622 else {
623 nlaydown[lay]++;
624 nlay[lay]++;
625 ndetdown[det]++;
626 ndet[det]++;
627 }
628 }
629
630 // checks minimum values
631 Bool_t isok=kTRUE;
632 for (Int_t j=0; j<6; j++) {
633 if (nlayup[j]<fNReqLayUp[j]) isok=kFALSE;
634 if (nlaydown[j]<fNReqLayDown[j]) isok=kFALSE;
635 if (nlay[j]<fNReqLay[j]) isok=kFALSE;
636 }
637 for (Int_t j=0; j<3; j++) {
638 if (ndetup[j]<fNReqDetUp[j]) isok=kFALSE;
639 if (ndetdown[j]<fNReqDetDown[j]) isok=kFALSE;
640 if (ndet[j]<fNReqDet[j]) isok=kFALSE;
641 }
642 if (!isok) {
75d480f6 643 AliDebug(2,Form("Track does not meet all location point requirements!"));
d3603b4d 644 return NULL;
645 }
646 }
647
648 // build a new track with (sorted) (prealigned) good points
649 atps=new AliTrackPointArray(ngoodpts);
650
651 for (int i=0; i<npts; i++) idx[i]=i;
652 // sort track if required
653 if (fUseSortedTracks) TMath::Sort(npts,atp->GetY(),idx); // sort descending...
654
655 Int_t npto=0;
656 for (int i=0; i<npts; i++) {
657 // skip not defined points
658 if (intidx[idx[i]]<0) continue;
659 atp->GetPoint(p,idx[i]);
660
661 // prealign point if required
662 // get IDEAL matrix
663 TGeoHMatrix *svOrigMatrix = fMilleModule[intidx[idx[i]]]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
664 // get back real local coordinates: use OriginalGlobalMatrix because AliTrackPoints were written
665 // with idel geometry
666 Double_t pg[3],pl[3];
667 pg[0]=p.GetX();
668 pg[1]=p.GetY();
669 pg[2]=p.GetZ();
670 AliDebug(3,Form("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
671 svOrigMatrix->MasterToLocal(pg,pl);
672
673 AliDebug(3,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]));
674
675 // update covariance matrix
676 TGeoHMatrix hcov;
677 Double_t hcovel[9];
678 hcovel[0]=double(p.GetCov()[0]);
679 hcovel[1]=double(p.GetCov()[1]);
680 hcovel[2]=double(p.GetCov()[2]);
681 hcovel[3]=double(p.GetCov()[1]);
682 hcovel[4]=double(p.GetCov()[3]);
683 hcovel[5]=double(p.GetCov()[4]);
684 hcovel[6]=double(p.GetCov()[2]);
685 hcovel[7]=double(p.GetCov()[4]);
686 hcovel[8]=double(p.GetCov()[5]);
687 hcov.SetRotation(hcovel);
688 // now rotate in local system
689 hcov.Multiply(svOrigMatrix);
690 hcov.MultiplyLeft(&svOrigMatrix->Inverse());
691 // now hcov is LOCAL COVARIANCE MATRIX
692
693
694 // pepopepo
695 if (fBug==1) {
696 // correzione bug LAYER 5 SSD temporanea..
697 int ssdidx=AliITSAlignMilleModule::GetIndexFromVolumeID(p.GetVolumeID());
698 if (ssdidx>=500 && ssdidx<1248) {
699 int ladder=(ssdidx-500)%22;
700 if (ladder==18) p.SetVolumeID(AliITSAlignMilleModule::GetVolumeIDFromIndex(ssdidx+1));
701 if (ladder==19) p.SetVolumeID(AliITSAlignMilleModule::GetVolumeIDFromIndex(ssdidx-1));
702 }
703 }
704
705 /// get (evenctually prealigned) matrix of sens. vol.
706 TGeoHMatrix *svMatrix = fMilleModule[intidx[idx[i]]]->GetSensitiveVolumeMatrix(p.GetVolumeID());
707 // modify global coordinates according with pre-aligment
708 svMatrix->LocalToMaster(pl,pg);
709 // now rotate in local system
710 hcov.Multiply(&svMatrix->Inverse());
711 hcov.MultiplyLeft(svMatrix);
712 // hcov is back in GLOBAL RF
713 Float_t pcov[6];
714 pcov[0]=hcov.GetRotationMatrix()[0];
715 pcov[1]=hcov.GetRotationMatrix()[1];
716 pcov[2]=hcov.GetRotationMatrix()[2];
717 pcov[3]=hcov.GetRotationMatrix()[4];
718 pcov[4]=hcov.GetRotationMatrix()[5];
719 pcov[5]=hcov.GetRotationMatrix()[8];
720
721 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
722 AliDebug(3,Form("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]));
723 atps->AddPoint(npto,&p);
724 AliDebug(2,Form("Adding point[%d] = ( %f , %f , %f ) volid = %d",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] ));
725
726 npto++;
727 }
728
729 return atps;
730}
731
732
733
734AliTrackPointArray *AliITSAlignMille::SortTrack(AliTrackPointArray *atp) {
735 /// sort alitrackpoints w.r.t. global Y direction
736 AliTrackPointArray *atps=NULL;
737 Int_t idx[20];
738 Int_t npts=atp->GetNPoints();
739 AliTrackPoint p;
740 atps=new AliTrackPointArray(npts);
741
742 TMath::Sort(npts,atp->GetY(),idx);
743
744 for (int i=0; i<npts; i++) {
745 atp->GetPoint(p,idx[i]);
746 atps->AddPoint(i,&p);
747 AliDebug(2,Form("Point[%d] = ( %f , %f , %f ) volid = %d",i,atps->GetX()[i],atps->GetY()[i],atps->GetZ()[i],atps->GetVolumeID()[i] ));
748 }
749 return atps;
750}
751
752
483fd7d8 753Int_t AliITSAlignMille::InitModuleParams() {
754 /// initialize geometry parameters for a given detector
755 /// for current cluster (fCluster)
756 /// fGlobalInitParam[] is set as:
757 /// [tx,ty,tz,psi,theta,phi]
483fd7d8 758 ///
483fd7d8 759 /// return 0 if success
760
761 if (!fGeoManager) {
762 AliInfo("ITS geometry not initialized!");
763 return -1;
764 }
765
0856765c 766 // now 'voluid' is the volumeID of a SENSITIVE VOLUME (coming from a cluster)
767
483fd7d8 768 // set the internal index (index in module list)
769 UShort_t voluid=fCluster.GetVolumeID();
770 Int_t k=fNModules-1;
d3603b4d 771 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
483fd7d8 772 if (k<0) return -3;
0856765c 773 fCurrentModuleInternalIndex=k; // the internal index of the SUPERMODULE
774
775 fCurrentModuleIndex=fMilleModule[k]->GetIndex(); // index of the SUPERMODULE
483fd7d8 776
777 // set the index
0856765c 778 Int_t index = GetModuleIndex(voluid);
483fd7d8 779 if (index<0) return -2;
0856765c 780 fCurrentSensVolIndex = index; // the index of the SENSITIVE VOLUME
483fd7d8 781
782 fModuleInitParam[0] = 0.0;
783 fModuleInitParam[1] = 0.0;
784 fModuleInitParam[2] = 0.0;
785 fModuleInitParam[3] = 0.0; // psi (X)
786 fModuleInitParam[4] = 0.0; // theta (Y)
787 fModuleInitParam[5] = 0.0; // phi (Z)
0856765c 788
0856765c 789 fCurrentModuleHMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetMatrix();
790
791 for (int ii=0; ii<3; ii++)
792 fCurrentModuleTranslation[ii]=fCurrentModuleHMatrix->GetTranslation()[ii];
793
d3603b4d 794 /// get (evenctually prealigned) matrix of sens. vol.
795 TGeoHMatrix *svMatrix = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeMatrix(voluid);
796
483fd7d8 797 fMeasGlo[0] = fCluster.GetX();
798 fMeasGlo[1] = fCluster.GetY();
d3603b4d 799 fMeasGlo[2] = fCluster.GetZ();
800 svMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
0856765c 801 AliDebug(2,Form("Local coordinates of measured point : X=%f Y=%f Z=%f \n",fMeasLoc[0] ,fMeasLoc[1] ,fMeasLoc[2] ));
0856765c 802
483fd7d8 803 // set stdev from cluster
804 TGeoHMatrix hcov;
805 Double_t hcovel[9];
806 hcovel[0]=double(fCluster.GetCov()[0]);
807 hcovel[1]=double(fCluster.GetCov()[1]);
d3603b4d 808 hcovel[2]=double(fCluster.GetCov()[2]);
483fd7d8 809 hcovel[3]=double(fCluster.GetCov()[1]);
d3603b4d 810 hcovel[4]=double(fCluster.GetCov()[3]);
483fd7d8 811 hcovel[5]=double(fCluster.GetCov()[4]);
d3603b4d 812 hcovel[6]=double(fCluster.GetCov()[2]);
483fd7d8 813 hcovel[7]=double(fCluster.GetCov()[4]);
814 hcovel[8]=double(fCluster.GetCov()[5]);
815 hcov.SetRotation(hcovel);
816 // now rotate in local system
0856765c 817 hcov.Multiply(svMatrix);
d3603b4d 818 hcov.MultiplyLeft(&svMatrix->Inverse());
483fd7d8 819
d3603b4d 820 // set local sigmas
483fd7d8 821 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
822 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4]));
823 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
824
d3603b4d 825 // set minimum value for SigmaLoc to 10 micron
0856765c 826 if (fSigmaLoc[0]<0.0010) fSigmaLoc[0]=0.0010;
827 if (fSigmaLoc[2]<0.0010) fSigmaLoc[2]=0.0010;
828
d3603b4d 829 // multiply local sigmas by global factor
830 fSigmaLoc[0] *= fSigmaXfactor;
831 fSigmaLoc[2] *= fSigmaZfactor;
832
833 // multiply local sigmas by individual factor
834 fSigmaLoc[0] *= fSensVolSigmaXfactor[index];
835 fSigmaLoc[2] *= fSensVolSigmaZfactor[index];
836
837 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%g fSigmaLocY=%g fSigmaLocZ=%g \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
483fd7d8 838
839 return 0;
840}
841
842void AliITSAlignMille::SetCurrentModule(Int_t index) {
0856765c 843 /// set as current the SuperModule that contains the 'index' sens.vol.
844 if (index<0 || index>2197) {
845 AliInfo("index does not correspond to a sensitive volume!");
846 return;
847 }
483fd7d8 848 UShort_t voluid=GetModuleVolumeID(index);
0856765c 849 Int_t k=IsContained(voluid);
850 if (k>=0){
0856765c 851 fCluster.SetVolumeID(voluid);
0856765c 852 fCluster.SetXYZ(0,0,0);
853 InitModuleParams();
854 }
855 else
d3603b4d 856 AliInfo(Form("module %d not defined\n",index));
0856765c 857}
858
859void AliITSAlignMille::SetCurrentSensitiveModule(Int_t index) {
860 /// set as current the SuperModule that contains the 'index' sens.vol.
861 if (index<0 || index>2197) {
862 AliInfo("index does not correspond to a sensitive volume!");
863 return;
864 }
865 UShort_t voluid=AliITSAlignMilleModule::GetVolumeIDFromIndex(index);
866 Int_t k=IsDefined(voluid);
0856765c 867 if (k>=0){
483fd7d8 868 fCluster.SetVolumeID(voluid);
869 fCluster.SetXYZ(0,0,0);
870 InitModuleParams();
871 }
0856765c 872 else
d3603b4d 873 AliInfo(Form("module %d not defined\n",index));
483fd7d8 874}
875
0856765c 876void AliITSAlignMille::Print(Option_t*) const
877{
d3603b4d 878 /// print infos
483fd7d8 879 printf("*** AliMillepede for ITS ***\n");
0856765c 880 printf(" number of defined super modules: %d\n",fNModules);
881
483fd7d8 882 if (fGeoManager)
883 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
884 else
885 printf(" geometry not loaded\n");
0856765c 886
887 if (fUseSuperModules)
888 printf(" using custom supermodules ( %d defined )\n",fNSuperModules);
889 else
890 printf(" custom supermodules not used\n");
891
892 if (fUsePreAlignment)
893 printf(" using prealignment from %s \n",fPreAlignmentFileName.Data());
894 else
895 printf(" prealignment not used\n");
896
d3603b4d 897 if (fBOn)
898 printf(" B Field set to %f T - using Riemann's helices\n",fBField);
899 else
900 printf(" B Field OFF - using straight lines \n");
901
483fd7d8 902 if (fUseLocalShifts)
903 printf(" Alignment shifts will be computed in LOCAL RS\n");
904 else
905 printf(" Alignment shifts will be computed in GLOBAL RS\n");
d3603b4d 906
907 if (fRequirePoints) printf(" Required points in tracks:\n");
908 for (Int_t i=0; i<6; i++) {
909 if (fNReqLayUp[i]>0) printf(" Layer %d : %d points with Y>0\n",i+1,fNReqLayUp[i]);
910 if (fNReqLayDown[i]>0) printf(" Layer %d : %d points with Y<0\n",i+1,fNReqLayDown[i]);
911 if (fNReqLay[i]>0) printf(" Layer %d : %d points \n",i+1,fNReqLay[i]);
912 }
913 for (Int_t i=0; i<3; i++) {
914 if (fNReqDetUp[i]>0) printf(" Detector %d : %d points with Y>0\n",i+1,fNReqDetUp[i]);
915 if (fNReqDetDown[i]>0) printf(" Detector %d : %d points with Y<0\n",i+1,fNReqDetDown[i]);
916 if (fNReqDet[i]>0) printf(" Detector %d : %d points \n",i+1,fNReqDet[i]);
917 }
0856765c 918
d3603b4d 919 printf("\n Millepede configuration parameters:\n");
920 printf(" init parsig for translations : %.4f\n",fParSigTranslations);
921 printf(" init parsig for rotations : %.4f\n",fParSigRotations);
922 printf(" init value for chi2 cut : %.4f\n",fStartFac);
923 printf(" first iteration cut value : %.4f\n",fResCutInitial);
924 printf(" other iterations cut value : %.4f\n",fResCut);
925 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
926 printf(" mult. fact. for local sigmas : %.4f %.4f\n",fSigmaXfactor, fSigmaZfactor);
483fd7d8 927
928 printf("List of defined modules:\n");
0856765c 929 printf(" intidx\tindex\tvoluid\tname\n");
483fd7d8 930 for (int i=0; i<fNModules; i++)
0856765c 931 printf(" %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],fMilleModule[i]->GetName());
932
483fd7d8 933}
934
0856765c 935AliITSAlignMilleModule *AliITSAlignMille::GetMilleModule(UShort_t voluid)
936{
937 // return pointer to a define supermodule
938 // return NULL if error
939 Int_t i=IsDefined(voluid);
940 if (i<0) return NULL;
941 return fMilleModule[i];
942}
943
75d480f6 944AliITSAlignMilleModule *AliITSAlignMille::GetCurrentModule() const
0856765c 945{
946 if (fNModules) return fMilleModule[fCurrentModuleInternalIndex];
947 return NULL;
948}
949
950void AliITSAlignMille::PrintCurrentModuleInfo()
951{
483fd7d8 952 ///
0856765c 953 Int_t k=fCurrentModuleInternalIndex;
954 if (k<0 || k>=fNModules) return;
955 fMilleModule[k]->Print("");
483fd7d8 956}
957
d3603b4d 958Bool_t AliITSAlignMille::InitRiemanFit() {
959 // Initialize Riemann Fitter for current track
960 // return kFALSE if error
961
962 if (!fBOn) return kFALSE;
963
964 Int_t npts=0;
965 AliTrackPoint ap;
966 npts = fTrack->GetNPoints();
967 AliDebug(3,Form("Fitting track with %d points",npts));
968
969 fRieman->Reset();
970 fRieman->SetTrackPointArray(fTrack);
971
972 TArrayI ai(npts);
973 for (Int_t ipt=0; ipt<npts; ipt++) ai[ipt]=fTrack->GetVolumeID()[ipt];
974
975 // fit track with 5 params in his own tracking-rotated reference system
976 // xc = -p[1]/p[0];
977 // yc = 1/p[0];
978 // R = sqrt( x0*x0 + y0*y0 - y0*p[2]);
979 if (!fRieman->Fit(&ai,NULL,(AliGeomManager::ELayerID)1,(AliGeomManager::ELayerID)6)) {
980 return kFALSE;
981 }
982
983 for (int i=0; i<5; i++)
984 fLocalInitParam[i] = fRieman->GetParam()[i];
985
986 return kTRUE;
987}
988
483fd7d8 989
990void AliITSAlignMille::InitTrackParams(int meth) {
991 /// initialize local parameters with different methods
992 /// for current track (fTrack)
993
0856765c 994 Int_t npts=0;
995 TF1 *f1=NULL;
996 TGraph *g=NULL;
997 Float_t sigmax[20],sigmay[20],sigmaz[20];
998 AliTrackPoint ap;
999 TGraphErrors *ge=NULL;
1000
483fd7d8 1001 switch (meth) {
1002 case 1: // simple linear interpolation
1003 // get local starting parameters (to be substituted by ESD track parms)
1004 // local parms (fLocalInitParam[]) are:
1005 // [0] = global x coord. of straight line intersection at y=0 plane
1006 // [1] = global z coord. of straight line intersection at y=0 plane
1007 // [2] = px/py
1008 // [3] = pz/py
1009
1010 // test #1: linear fit in x(y) and z(y)
0856765c 1011 npts = fTrack->GetNPoints();
d3603b4d 1012 AliDebug(3,Form("*** initializing track with %d points ***",npts));
483fd7d8 1013
0856765c 1014 f1=new TF1("f1","[0]+x*[1]",-50,50);
483fd7d8 1015
0856765c 1016 g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
483fd7d8 1017 g->Fit(f1,"RNQ");
1018 fLocalInitParam[0] = f1->GetParameter(0);
1019 fLocalInitParam[2] = f1->GetParameter(1);
1020 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
1021 delete g; g=NULL;
1022
1023 g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ());
1024 g->Fit(f1,"RNQ");
1025 fLocalInitParam[1] = f1->GetParameter(0);
1026 fLocalInitParam[3] = f1->GetParameter(1);
1027 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
1028 delete g;
1029 delete f1;
1030
1031 break;
0856765c 1032
1033 case 2: // simple linear interpolation weighted using sigmas
1034 // get local starting parameters (to be substituted by ESD track parms)
1035 // local parms (fLocalInitParam[]) are:
1036 // [0] = global x coord. of straight line intersection at y=0 plane
1037 // [1] = global z coord. of straight line intersection at y=0 plane
1038 // [2] = px/py
1039 // [3] = pz/py
1040
1041 // test #1: linear fit in x(y) and z(y)
1042 npts = fTrack->GetNPoints();
d3603b4d 1043 AliDebug(3,Form("*** initializing track with %d points ***",npts));
1044
0856765c 1045 for (Int_t isig=0; isig<npts; isig++) {
1046 fTrack->GetPoint(ap,isig);
1047 sigmax[isig]=ap.GetCov()[0];
1048 if (sigmax[isig]<1.0e-07) sigmax[isig]=1.0e-07; // minimum sigma=3 mu
1049 sigmax[isig]=TMath::Sqrt(sigmax[isig]);
1050
d3603b4d 1051 sigmay[isig]=ap.GetCov()[3];
0856765c 1052 if (sigmay[isig]<1.0e-07) sigmay[isig]=1.0e-07; // minimum sigma=3 mu
1053 sigmay[isig]=TMath::Sqrt(sigmay[isig]);
1054
1055 sigmaz[isig]=ap.GetCov()[5];
1056 if (sigmaz[isig]<1.0e-07) sigmaz[isig]=1.0e-07; // minimum sigma=3 mu
1057 sigmaz[isig]=TMath::Sqrt(sigmaz[isig]);
d3603b4d 1058
1059 if (fTempExcludedModule>=0 && fTempExcludedModule==AliITSAlignMilleModule::GetIndexFromVolumeID(ap.GetVolumeID())) {
1060 sigmax[isig] *= 1000.;
1061 sigmay[isig] *= 1000.;
1062 sigmaz[isig] *= 1000.;
1063 fTempExcludedModule=-1;
1064 }
0856765c 1065 }
1066
1067 f1=new TF1("f1","[0]+x*[1]",-50,50);
1068
1069 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetX(),sigmay,sigmax);
1070 ge->Fit(f1,"RNQ");
1071 fLocalInitParam[0] = f1->GetParameter(0);
1072 fLocalInitParam[2] = f1->GetParameter(1);
1073 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
1074 delete ge; ge=NULL;
1075
1076 ge=new TGraphErrors(npts,fTrack->GetY(),fTrack->GetZ(),sigmay,sigmaz);
1077 ge->Fit(f1,"RNQ");
1078 fLocalInitParam[1] = f1->GetParameter(0);
1079 fLocalInitParam[3] = f1->GetParameter(1);
1080 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
1081 delete ge;
1082 delete f1;
1083
1084 break;
1085
483fd7d8 1086 }
0856765c 1087}
483fd7d8 1088
0856765c 1089Int_t AliITSAlignMille::IsDefined(UShort_t voluid) const
1090{
1091 // checks if supermodule 'voluid' is defined and return the internal index
1092 // return -1 if error
1093 Int_t k=fNModules-1;
1094 while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
1095 if (k<0) return -1;
1096 return k;
1097}
1098
1099Int_t AliITSAlignMille::IsContained(UShort_t voluid) const
1100{
1101 // checks if the sensitive module 'voluid' is contained inside a supermodule and return the internal index of the last identified supermodule
1102 // return -1 if error
1103 if (AliITSAlignMilleModule::GetIndexFromVolumeID(voluid)<0) return -1;
1104 Int_t k=fNModules-1;
1105 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
1106 if (k<0) return -1;
1107 return k;
483fd7d8 1108}
0856765c 1109
483fd7d8 1110Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
1111{
0856765c 1112 /// check if a sensitive volume is contained inside one of the defined supermodules
483fd7d8 1113 Int_t k=fNModules-1;
0856765c 1114 while (k>=0 && !(fMilleModule[k]->IsIn(voluid)) ) k--;
483fd7d8 1115 if (k>=0) return kTRUE;
1116 return kFALSE;
1117}
1118
1119Int_t AliITSAlignMille::CheckCurrentTrack() {
1120 /// checks if AliTrackPoints belongs to defined modules
1121 /// return number of good poins
1122 /// return 0 if not enough points
1123
1124 Int_t npts = fTrack->GetNPoints();
1125 Int_t ngoodpts=0;
1126 // debug points
1127 for (int j=0; j<npts; j++) {
1128 UShort_t voluid = fTrack->GetVolumeID()[j];
1129 if (CheckVolumeID(voluid)) {
1130 ngoodpts++;
1131 }
1132 }
d3603b4d 1133
483fd7d8 1134 if (ngoodpts<fMinNPtsPerTrack) return 0;
1135
1136 return ngoodpts;
1137}
1138
1139Int_t AliITSAlignMille::ProcessTrack(AliTrackPointArray *track) {
1140 /// Process track; Loop over hits and set local equations
1141 /// here 'track' is a AliTrackPointArray
1142 /// return 0 if success;
1143
1144 if (!fIsMilleInit) {
1145 AliInfo("Millepede has not been initialized!");
1146 return -1;
1147 }
1148
1149 Int_t npts = track->GetNPoints();
d3603b4d 1150 AliDebug(2,Form("*** Input track with %d points ***",npts));
483fd7d8 1151
d3603b4d 1152 // preprocessing of the input track: keep only points in defined volumes,
1153 // move points if prealignment is set, sort by Yglo if required
1154
1155 fTrack=PrepareTrack(track);
1156 if (!fTrack) return -1;
483fd7d8 1157
d3603b4d 1158 npts = fTrack->GetNPoints();
1159 AliDebug(2,Form("*** Processing prepared track with %d points ***",npts));
1160
1161 if (!fBOn) { // straight lines
1162 // set local starting parameters (to be substituted by ESD track parms)
1163 // local parms (fLocalInitParam[]) are:
1164 // [0] = global x coord. of straight line intersection at y=0 plane
1165 // [1] = global z coord. of straight line intersection at y=0 plane
1166 // [2] = px/py
1167 // [3] = pz/py
1168 InitTrackParams(fInitTrackParamsMeth);
1169 }
1170 else {
1171 // local parms (fLocalInitParam[]) are the Riemann Fitter params
1172 if (!InitRiemanFit()) {
1173 AliInfo("Riemann fit failed! skipping this track...");
1174 delete fTrack;
1175 fTrack=NULL;
1176 return -5;
1177 }
1178 }
1179
1180 Int_t nloceq=0;
75d480f6 1181 AliITSAlignMilleData *md = new AliITSAlignMilleData[npts];
d3603b4d 1182
483fd7d8 1183 for (Int_t ipt=0; ipt<npts; ipt++) {
1184 fTrack->GetPoint(fCluster,ipt);
d3603b4d 1185 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
483fd7d8 1186
1187 // set geometry parameters for the the current module
483fd7d8 1188 if (InitModuleParams()) continue;
483fd7d8 1189 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
d3603b4d 1190 AliDebug(2,Form(" Preprocessed Point = ( %f , %f , %f ) \n",fCluster.GetX(),fCluster.GetY(),fCluster.GetZ()));
1191
1192 if (!AddLocalEquation(md[nloceq])) {
1193 nloceq++;
1194 fProcessedPoints[fCurrentModuleInternalIndex]++;
1195 }
1196 else {
1197 fTotBadLocEqPoints++;
1198 }
483fd7d8 1199
483fd7d8 1200 } // end loop over points
1201
d3603b4d 1202 delete fTrack;
1203 fTrack=NULL;
1204
1205 // not enough good points!
1206 if (nloceq<fMinNPtsPerTrack) {
1207 delete [] md;
1208 return -1;
1209 }
1210
1211 // finally send local equations to millepede
1212 SetLocalEquations(md,nloceq);
1213
1214 delete [] md;
1215
483fd7d8 1216 return 0;
1217}
1218
1219Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) {
d3603b4d 1220 /// calculate intersection point of track with current module in local coordinates
1221 /// according with a given set of parameters (local(4/5) and global(6))
483fd7d8 1222 /// and fill fPintLoc/Glo
d3603b4d 1223 /// local are: pgx0, pgz0, ugx, ugz OR riemann fitters pars
483fd7d8 1224 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
1225 /// return 0 if success
1226
d3603b4d 1227 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]));
1228 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 1229
483fd7d8 1230
d3603b4d 1231 // prepare the TGeoHMatrix
0856765c 1232 TGeoHMatrix *fTempHMat = fMilleModule[fCurrentModuleInternalIndex]->GetSensitiveVolumeModifiedMatrix(fCluster.GetVolumeID(),gpar);
1233 if (!fTempHMat) return -1;
d3603b4d 1234
1235 Double_t v0g[3]; // vector with straight line direction in global coord.
1236 Double_t p0g[3]; // point of the straight line (glo)
1237
1238 if (fBOn) { // B FIELD!
1239 AliTrackPoint prf;
1240 for (int ip=0; ip<5; ip++)
1241 fRieman->SetParam(ip,lpar[ip]);
1242
1243 if (!fRieman->GetPCA(fCluster,prf)) {
1244 AliInfo(Form("error in GetPCA for point %d",fCluster.GetVolumeID()));
1245 return -3;
1246 }
1247 // now determine straight line passing tangent to fit curve at prf
1248 // ugx = dX/dY_glo = DeltaX/DeltaY_glo
1249 // mo' P1=(X,Y,Z)_glo_prf
1250 // => (x,y,Z)_trk_prf ruotando di alpha...
1251 Double_t alpha=fRieman->GetAlpha();
1252 Double_t x1g = prf.GetX();
1253 Double_t y1g = prf.GetY();
1254 Double_t z1g = prf.GetZ();
1255 Double_t x1t = x1g*TMath::Cos(alpha) + y1g*TMath::Sin(alpha);
1256 Double_t y1t = -x1g*TMath::Sin(alpha) + y1g*TMath::Cos(alpha);
1257 Double_t z1t = z1g;
1258
1259 Double_t x2t = x1t+1.0;
1260 Double_t y2t = y1t+fRieman->GetDYat(x1t);
1261 Double_t z2t = z1t+fRieman->GetDZat(x1t);
1262 Double_t x2g = x2t*TMath::Cos(alpha) - y2t*TMath::Sin(alpha);
1263 Double_t y2g = x2t*TMath::Sin(alpha) + y2t*TMath::Cos(alpha);
1264 Double_t z2g = z2t;
1265
1266 AliDebug(3,Form("Riemann frame: fAlpha = %f = %f ",alpha,alpha*180./TMath::Pi()));
1267 AliDebug(3,Form(" prf_glo=( %f , %f , %f ) prf_rf=( %f , %f , %f )\n", x1g,y1g,z1g, x1t,y1t,z1t));
1268 AliDebug(3,Form(" mov_glo=( %f , %f , %f ) rf=( %f , %f , %f )\n",x2g,y2g,z2g, x2t,y2t,z2t));
1269
1270 if (TMath::Abs(y2g-y1g)<1e-15) {
1271 AliInfo("DeltaY=0! Cannot proceed...");
1272 return -1;
1273 }
1274 // ugx,1,ugz
1275 v0g[0] = (x2g-x1g)/(y2g-y1g);
1276 v0g[1]=1.0;
1277 v0g[2] = (z2g-z1g)/(y2g-y1g);
1278
1279 // point: just keep prf
1280 p0g[0]=x1g;
1281 p0g[1]=y1g;
1282 p0g[2]=z1g;
1283 }
1284 else { // staight line
1285 // vector of initial straight line direction in glob. coord
1286 v0g[0]=lpar[2];
1287 v0g[1]=1.0;
1288 v0g[2]=lpar[3];
1289
1290 // intercept in yg=0 plane in glob coord
1291 p0g[0]=lpar[0];
1292 p0g[1]=0.0;
1293 p0g[2]=lpar[1];
1294 }
1295 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]));
1296
483fd7d8 1297 // same in local coord.
1298 Double_t p0l[3],v0l[3];
1299 fTempHMat->MasterToLocalVect(v0g,v0l);
1300 fTempHMat->MasterToLocal(p0g,p0l);
d3603b4d 1301
483fd7d8 1302 if (TMath::Abs(v0l[1])<1e-15) {
1303 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
1304 return -1;
1305 }
1306
1307 // local intersection point
1308 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
1309 fPintLoc[1] = 0;
1310 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
1311
1312 // global intersection point
1313 fTempHMat->LocalToMaster(fPintLoc,fPintGlo);
1314 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 1315
483fd7d8 1316 return 0;
1317}
1318
1319Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) {
1320 /// calculate numerically (ROOT's style) the derivatives for
1321 /// local X intersection and local Z intersection
d3603b4d 1322 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0 OR riemann's params
483fd7d8 1323 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
1324 /// return 0 if success
1325
1326 // copy initial parameters
75d480f6 1327 Double_t lpar[ITSMILLENLOCAL];
1328 Double_t gpar[ITSMILLENPARCH];
1329 for (Int_t i=0; i<ITSMILLENLOCAL; i++) lpar[i]=fLocalInitParam[i];
1330 for (Int_t i=0; i<ITSMILLENPARCH; i++) gpar[i]=fModuleInitParam[i];
483fd7d8 1331
483fd7d8 1332 // trial with fixed dpar...
1333 Double_t dpar=0.0;
d3603b4d 1334
1335 if (islpar) { // track parameters
483fd7d8 1336 //dpar=fLocalInitParam[paridx]*0.001;
1337 // set minimum dpar
d3603b4d 1338 if (!fBOn) {
1339 if (paridx<2) dpar=1.0e-4; // translations
1340 else dpar=1.0e-6; // direction
1341 }
1342 else { // B Field
1343 // pepo: proviamo con 1/1000, poi evenctually 1/100...
1344 Double_t dfrac=0.01;
1345 switch(paridx) {
1346 case 0:
1347 // RMS cosmics: 1e-4
1348 dpar = TMath::Max(1.0e-6,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1349 break;
1350 case 1:
1351 // RMS cosmics: 0.2
1352 dpar = TMath::Max(0.002,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1353 break;
1354 case 2:
1355 // RMS cosmics: 9
1356 dpar = TMath::Max(0.09,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1357 break;
1358 case 3:
1359 // RMS cosmics: 7
1360 dpar = TMath::Max(0.07,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1361 break;
1362 case 4:
1363 // RMS cosmics: 0.3
1364 dpar = TMath::Max(0.003,TMath::Abs(fLocalInitParam[paridx]*dfrac));
1365 break;
1366 }
1367 }
483fd7d8 1368 }
d3603b4d 1369 else { // alignment global parameters
483fd7d8 1370 //dpar=fModuleInitParam[paridx]*0.001;
1371 if (paridx<3) dpar=1.0e-4; // translations
1372 else dpar=1.0e-2; // angles
1373 }
d3603b4d 1374
1375 AliDebug(3,Form("+++ using dpar=%g",dpar));
483fd7d8 1376
1377 // calculate derivative ROOT's like:
1378 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
1379 Double_t pintl1[3]; // f(x-h)
1380 Double_t pintl2[3]; // f(x-h/2)
1381 Double_t pintl3[3]; // f(x+h/2)
1382 Double_t pintl4[3]; // f(x+h)
1383
1384 // first values
1385 if (islpar) lpar[paridx] -= dpar;
1386 else gpar[paridx] -= dpar;
1387 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1388 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
1389
1390 // second values
1391 if (islpar) lpar[paridx] += dpar/2;
1392 else gpar[paridx] += dpar/2;
1393 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1394 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
1395
1396 // third values
1397 if (islpar) lpar[paridx] += dpar;
1398 else gpar[paridx] += dpar;
1399 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1400 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
1401
1402 // fourth values
1403 if (islpar) lpar[paridx] += dpar/2;
1404 else gpar[paridx] += dpar/2;
1405 if (CalcIntersectionPoint(lpar, gpar)) return -2;
1406 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
1407
1408 Double_t h2 = 1./(2.*dpar);
1409 Double_t d0 = pintl4[0]-pintl1[0];
1410 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
1411 fDerivativeXLoc = h2*(4*d2 - d0)/3.;
1412 if (TMath::Abs(fDerivativeXLoc) < 1.0e-9) fDerivativeXLoc=0.0;
1413
1414 d0 = pintl4[2]-pintl1[2];
1415 d2 = 2.*(pintl3[2]-pintl2[2]);
1416 fDerivativeZLoc = h2*(4*d2 - d0)/3.;
1417 if (TMath::Abs(fDerivativeZLoc) < 1.0e-9) fDerivativeZLoc=0.0;
1418
1419 AliDebug(3,Form("\n+++ derivatives +++ \n"));
1420 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",fDerivativeXLoc));
1421 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",fDerivativeZLoc));
1422
1423 return 0;
1424}
1425
d3603b4d 1426
75d480f6 1427Int_t AliITSAlignMille::AddLocalEquation(AliITSAlignMilleData &m) {
d3603b4d 1428 /// Define local equation for current cluster in X and Z coor.
1429 /// and store them to memory
483fd7d8 1430 /// return 0 if success
1431
1432 // store first interaction point
d3603b4d 1433 if (CalcIntersectionPoint(fLocalInitParam, fModuleInitParam)) return -4;
483fd7d8 1434 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
0856765c 1435 AliDebug(2,Form("Intesect. point: L( %f , %f , %f )",fPintLoc[0],fPintLoc[1],fPintLoc[2]));
483fd7d8 1436
1437 // calculate local derivatives numerically
75d480f6 1438 Double_t dXdL[ITSMILLENLOCAL],dZdL[ITSMILLENLOCAL];
d3603b4d 1439 for (Int_t i=0; i<fNLocal; i++) {
483fd7d8 1440 if (CalcDerivatives(i,kTRUE)) return -1;
1441 dXdL[i]=fDerivativeXLoc;
1442 dZdL[i]=fDerivativeZLoc;
1443 }
1444
75d480f6 1445 Double_t dXdG[ITSMILLENPARCH],dZdG[ITSMILLENPARCH];
1446 for (Int_t i=0; i<ITSMILLENPARCH; i++) {
483fd7d8 1447 if (CalcDerivatives(i,kFALSE)) return -1;
1448 dXdG[i]=fDerivativeXLoc;
1449 dZdG[i]=fDerivativeZLoc;
1450 }
1451
1452 AliDebug(2,Form("\n***************\n"));
d3603b4d 1453 for (Int_t i=0; i<fNLocal; i++)
483fd7d8 1454 AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
75d480f6 1455 for (Int_t i=0; i<ITSMILLENPARCH; i++)
483fd7d8 1456 AliDebug(2,Form("Global parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdG[i],dZdG[i]));
1457 AliDebug(2,Form("\n***************\n"));
1458
d3603b4d 1459 // check if at least 1 local and 1 global derivs. are not null
1460 Double_t nonzero=0.0;
1461 for (Int_t i=0; i<fNLocal; i++) nonzero += TMath::Abs(dXdL[i]);
1462 if (nonzero==0.0) {
1463 AliInfo("Aborting local equations for this point beacuse of zero local X derivatives!");
1464 return -2;
1465 }
1466 nonzero=0.0;
75d480f6 1467 for (Int_t i=0; i<ITSMILLENPARCH; i++) nonzero += TMath::Abs(dXdG[i]);
d3603b4d 1468 if (nonzero==0.0) {
1469 AliInfo("Aborting local equations for this point beacuse of zero global X derivatives!");
1470 return -2;
1471 }
1472 nonzero=0.0;
1473 for (Int_t i=0; i<fNLocal; i++) nonzero += TMath::Abs(dZdL[i]);
1474 if (nonzero==0.0) {
1475 AliInfo("Aborting local equations for this point beacuse of zero local Z derivatives!");
1476 return -2;
1477 }
1478 nonzero=0.0;
75d480f6 1479 for (Int_t i=0; i<ITSMILLENPARCH; i++) nonzero += TMath::Abs(dZdG[i]);
d3603b4d 1480 if (nonzero==0.0) {
1481 AliInfo("Aborting local equations for this point beacuse of zero global Z derivatives!");
1482 return -2;
1483 }
1484
1485 // ok, can copy to m
1486
1487 AliDebug(2,Form("Adding local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
483fd7d8 1488 // set equation for Xloc coordinate
d3603b4d 1489 for (Int_t i=0; i<fNLocal; i++) {
75d480f6 1490 m.GetIdxlocX()[i]=i;
1491 m.GetDerlocX()[i]=dXdL[i];
d3603b4d 1492 }
75d480f6 1493 for (Int_t i=0; i<ITSMILLENPARCH; i++) {
1494 m.GetIdxgloX()[i]=fCurrentModuleInternalIndex*ITSMILLENPARCH+i;
1495 m.GetDergloX()[i]=dXdG[i];
d3603b4d 1496 }
75d480f6 1497 m.SetMeasX(fMeasLoc[0]-fPintLoc0[0]);
1498 m.SetSigmaX(fSigmaLoc[0]);
483fd7d8 1499
d3603b4d 1500 AliDebug(2,Form("Adding local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
483fd7d8 1501 // set equation for Zloc coordinate
d3603b4d 1502 for (Int_t i=0; i<fNLocal; i++) {
75d480f6 1503 m.GetIdxlocZ()[i]=i;
1504 m.GetDerlocZ()[i]=dZdL[i];
d3603b4d 1505 }
75d480f6 1506 for (Int_t i=0; i<ITSMILLENPARCH; i++) {
1507 m.GetIdxgloZ()[i]=fCurrentModuleInternalIndex*ITSMILLENPARCH+i;
1508 m.GetDergloZ()[i]=dZdG[i];
d3603b4d 1509 }
75d480f6 1510 m.SetMeasZ(fMeasLoc[2]-fPintLoc0[2]);
1511 m.SetSigmaZ(fSigmaLoc[2]);
d3603b4d 1512
483fd7d8 1513 return 0;
1514}
1515
1516
75d480f6 1517void AliITSAlignMille::SetLocalEquations(AliITSAlignMilleData *m, Int_t neq) {
d3603b4d 1518 /// Set local equations with data stored in m
1519 /// return 0 if success
1520
1521 for (Int_t j=0; j<neq; j++) {
1522
75d480f6 1523 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",m[j].GetMeasX(), m[j].GetSigmaX()));
d3603b4d 1524 // set equation for Xloc coordinate
1525 for (Int_t i=0; i<fNLocal; i++)
75d480f6 1526 SetLocalDerivative( m[j].GetIdxlocX()[i], m[j].GetDerlocX()[i] );
d3603b4d 1527
75d480f6 1528 for (Int_t i=0; i<ITSMILLENPARCH; i++)
1529 SetGlobalDerivative( m[j].GetIdxgloX()[i], m[j].GetDergloX()[i] );
d3603b4d 1530
75d480f6 1531 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m[j].GetMeasX(), m[j].GetSigmaX());
d3603b4d 1532
1533
75d480f6 1534 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",m[j].GetMeasZ(), m[j].GetSigmaZ()));
d3603b4d 1535 // set equation for Zloc coordinate
1536 for (Int_t i=0; i<fNLocal; i++)
75d480f6 1537 SetLocalDerivative( m[j].GetIdxlocZ()[i], m[j].GetDerlocZ()[i] );
d3603b4d 1538
75d480f6 1539 for (Int_t i=0; i<ITSMILLENPARCH; i++)
1540 SetGlobalDerivative( m[j].GetIdxgloZ()[i], m[j].GetDergloZ()[i] );
d3603b4d 1541
75d480f6 1542 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, m[j].GetMeasZ(), m[j].GetSigmaZ());
d3603b4d 1543 }
1544}
1545
1546
483fd7d8 1547void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
1548 /// Call local fit for this track
1549 if (!fIsMilleInit) {
1550 AliInfo("Millepede has not been initialized!");
1551 return;
1552 }
1553 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
1554 AliDebug(2,Form("iRes = %d",iRes));
75d480f6 1555 //if (iRes && !lSingleFit) {
1556 if (!lSingleFit) { // Ruben Shahoyan's bug fix
483fd7d8 1557 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
1558 }
1559}
1560
1561void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
1562 /// Call global fit; Global parameters are stored in parameters
1563 if (!fIsMilleInit) {
1564 AliInfo("Millepede has not been initialized!");
1565 return;
1566 }
1567 fMillepede->GlobalFit(parameters,errors,pulls);
1568 AliInfo("Done fitting global parameters!");
1569}
1570
1571Double_t AliITSAlignMille::GetParError(Int_t iPar) {
1572 /// Get error of parameter iPar
1573 if (!fIsMilleInit) {
1574 AliInfo("Millepede has not been initialized!");
1575 return 0;
1576 }
1577 Double_t lErr = fMillepede->GetParError(iPar);
1578 return lErr;
1579}
1580
1581void AliITSAlignMille::PrintGlobalParameters() {
1582 /// Print global parameters
1583 if (!fIsMilleInit) {
1584 AliInfo("Millepede has not been initialized!");
1585 return;
1586 }
1587 fMillepede->PrintGlobalParameters();
1588}
1589
1590// //_________________________________________________________________________
0856765c 1591Int_t AliITSAlignMille::LoadSuperModuleFile(const Char_t *sfile)
1592{
1593 // load definitions of supermodules from a root file
1594 // return 0 if success
1595
1596 TFile *smf=TFile::Open(sfile);
1597 if (!smf->IsOpen()) {
1598 AliInfo(Form("Cannot open supermodule file %s",sfile));
1599 return -1;
1600 }
1601
1602 TClonesArray *sma=(TClonesArray*)smf->Get("ITSMilleSuperModules");
1603 if (!sma) {
1604 AliInfo(Form("Cannot find ITSMilleSuperModules array in file"));
1605 return -2;
1606 }
1607 Int_t nsma=sma->GetEntriesFast();
1608 AliInfo(Form("Array of SuperModules with %d entries\n",nsma));
1609
41815654 1610 Char_t st[2048];
1611 char symname[250];
0856765c 1612 UShort_t volid;
1613 TGeoHMatrix m;
1614
1615 for (Int_t i=0; i<nsma; i++) {
1616 AliAlignObjParams *a = (AliAlignObjParams*)sma->UncheckedAt(i);
1617 volid=a->GetVolUID();
1618 strcpy(st,a->GetSymName());
1619 a->GetMatrix(m);
1620
1621 sscanf(st,"%s",symname);
1622 // decode module list
1623 char *stp=strstr(st,"ModuleList:");
1624 if (!stp) return -3;
1625 stp += 11;
1626 int idx[2200];
1627 char spp[200]; int jp=0;
1628 char cl[20];
1629 strcpy(st,stp);
1630 int l=strlen(st);
1631 int j=0;
1632 int n=0;
1633
1634 while (j<=l) {
1635 if (st[j]==9 || st[j]==32 || st[j]==10 || st[j]==0) {
1636 spp[jp]=0;
1637 jp=0;
1638 if (strlen(spp)) {
1639 int k=strcspn(spp,"-");
1640 if (k<int(strlen(spp))) { // c'e' il -
1641 strcpy(cl,&(spp[k+1]));
1642 spp[k]=0;
1643 int ifrom=atoi(spp); int ito=atoi(cl);
1644 for (int b=ifrom; b<=ito; b++) {
1645 idx[n]=b;
1646 n++;
1647 }
1648 }
1649 else { // numerillo singolo
1650 idx[n]=atoi(spp);
1651 n++;
1652 }
1653 }
1654 }
1655 else {
1656 spp[jp]=st[j];
1657 jp++;
1658 }
1659 j++;
1660 }
1661 UShort_t volidsv[2198];
1662 for (j=0;j<n;j++) {
1663 volidsv[j]=AliITSAlignMilleModule::GetVolumeIDFromIndex(idx[j]);
1664 if (!volidsv[j]) {
1665 AliInfo(Form("Index %d not valid (range 0->2197)",idx[j]));
1666 return -5;
1667 }
1668 }
1669 Int_t smindex=int(2198+volid-14336); // virtual index
1670 fSuperModule[fNSuperModules]=new AliITSAlignMilleModule(smindex,volid,symname,&m,n,volidsv);
1671
1672 //-------------
1673 fNSuperModules++;
1674 }
1675
1676 smf->Close();
1677
1678 fUseSuperModules=1;
1679 return 0;
1680}
483fd7d8 1681