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