]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITSAlignMille.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
New class for ITS alignment with Millepede (M. Lunardon)
[u/mrichter/AliRoot.git] / ITS / AliITSAlignMille.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18//-----------------------------------------------------------------------------
19/// \class AliITSAlignMille
20/// Alignment class fro the ALICE ITS detector
21///
22/// ITS specific alignment class which interface to AliMillepede.
23/// For each track ProcessTrack calculates the local and global derivatives
24/// at each hit and fill the corresponding local equations. Provide methods for
25/// fixing or constraining detection elements for best results.
26///
27/// \author M. Lunardon (thanks to J. Castillo)
28//-----------------------------------------------------------------------------
29
30#include <TF1.h>
31#include <TGraph.h>
32#include <TGeoMatrix.h>
33#include <TMath.h>
34
35#include "AliITSAlignMille.h"
36#include "AliITSgeomTGeo.h"
37#include "AliGeomManager.h"
38#include "AliMillepede.h"
39#include "AliTrackPointArray.h"
40#include "AliAlignObjParams.h"
41#include "AliLog.h"
42#include "TSystem.h" // come si fa?
43
44/// \cond CLASSIMP
45ClassImp(AliITSAlignMille)
46/// \endcond
47
48Int_t AliITSAlignMille::fgNDetElem = ITSMILLE_NDETELEM;
49Int_t AliITSAlignMille::fgNParCh = ITSMILLE_NPARCH;
50
51AliITSAlignMille::AliITSAlignMille(const Char_t *configFilename, Bool_t initmille)
52 : TObject(),
53 fMillepede(0),
54 fStartFac(16.),
55 fResCutInitial(100.),
56 fResCut(100.),
57 fNGlobal(ITSMILLE_NDETELEM*ITSMILLE_NPARCH),
58 fNLocal(ITSMILLE_NLOCAL),
59 fNStdDev(ITSMILLE_NSTDEV),
60 fIsMilleInit(kFALSE),
61 fParSigTranslations(0.0100),
62 fParSigRotations(0.1),
63 fTrack(NULL),
64 fCluster(),
65 fGlobalDerivatives(NULL),
66 fTempHMat(NULL),
67 fTempAlignObj(NULL),
68 fDerivativeXLoc(0),
69 fDerivativeZLoc(0),
70 fDeltaPar(0),
71 fMinNPtsPerTrack(3),
72 fGeometryFileName("geometry.root"),
73 fGeoManager(0),
74 fCurrentModuleIndex(0),
75 fCurrentModuleInternalIndex(0),
76 fNModules(0),
77 fUseLocalShifts(kTRUE),
78 fCurrentModuleHMatrix(NULL)
79{
80 /// main constructor that takes input from configuration file
81
82 fMillepede = new AliMillepede();
83 fGlobalDerivatives = new Double_t[fNGlobal];
84 fTempHMat = new TGeoHMatrix;
85 fCurrentModuleHMatrix = new TGeoHMatrix;
86
87 Int_t lc=LoadConfig(configFilename);
88 if (lc) {
89 AliInfo(Form("Error %d loading configuration from %s",lc,configFilename));
90 }
91 else {
92 if (initmille && fNGlobal<20000) {
93 AliInfo(Form("Initializing Millepede with %d gpar, %d lpar and %d stddev ...",fNGlobal, fNLocal, fNStdDev));
94 Init(fNGlobal, fNLocal, fNStdDev);
95 ResetLocalEquation();
96 AliInfo("Parameters initialized to zero");
97 }
98 else {
99 AliInfo("Millepede has not been initialized ... ");
100 }
101 }
102
103 fDeltaPar=0.0; // not used at the moment - to be checked later...
104
105}
106
107AliITSAlignMille::~AliITSAlignMille() {
108 /// Destructor
109 if (fMillepede) delete fMillepede;
110 delete [] fGlobalDerivatives;
111 delete fCurrentModuleHMatrix;
112 delete fTempHMat;
113}
114
115///////////////////////////////////////////////////////////////////////
116
117Int_t AliITSAlignMille::LoadConfig(const Char_t *cfile) {
118 /// return 0 if success
119 /// 1 if error in module index or voluid
120
121 FILE *pfc=fopen(cfile,"r");
122 if (!pfc) return -1;
123
124 Char_t st[200],st2[200];
125 Char_t tmp[100];
126 Int_t idx,itx,ity,itz,ith,ips,iph;
127 Float_t f1;
128 UShort_t voluid;
129 Int_t nmod=0;
130
131 while (fgets(st,200,pfc)) {
132
133 // skip comments
134 for (int i=0; i<int(strlen(st)); i++) {
135 if (st[i]=='#') st[i]=0;
136 }
137
138 if (strstr(st,"GEOMETRY_FILE")) {
139 sscanf(st,"%s %s",tmp,st2);
140 if (gSystem->AccessPathName(st2)) {
141 AliInfo("*** WARNING! *** geometry file not found! ");
142 return -1;
143 }
144 fGeometryFileName=st2;
145 InitGeometry();
146 }
147
148 if (strstr(st,"SET_PARSIG_TRA")) {
149 sscanf(st,"%s %f",tmp,&f1);
150 fParSigTranslations=f1;
151 }
152
153 if (strstr(st,"SET_PARSIG_ROT")) {
154 sscanf(st,"%s %f",tmp,&f1);
155 fParSigRotations=f1;
156 }
157
158 if (strstr(st,"SET_NSTDDEV")) {
159 sscanf(st,"%s %d",tmp,&idx);
160 fNStdDev=idx;
161 }
162
163 if (strstr(st,"SET_RESCUT_INIT")) {
164 sscanf(st,"%s %f",tmp,&f1);
165 fResCutInitial=f1;
166 }
167
168 if (strstr(st,"SET_RESCUT_OTHER")) {
169 sscanf(st,"%s %f",tmp,&f1);
170 fResCut=f1;
171 }
172
173 if (strstr(st,"SET_STARTFAC")) {
174 sscanf(st,"%s %f",tmp,&f1);
175 fStartFac=f1;
176 }
177
178 if (strstr(st,"SET_LOCAL_SHIFTS")) { // only enabled mode...
179 fUseLocalShifts = kTRUE;
180 }
181
182 if (strstr(st,"MODULE_INDEX")) {
183 sscanf(st,"%s %d %d %d %d %d %d %d",tmp,&idx,&itx,&ity,&itz,&iph,&ith,&ips);
184 voluid=GetModuleVolumeID(idx);
185 if (!voluid) return 1; // bad index
186 fModuleIndex[nmod]=idx;
187 fModuleVolumeID[nmod]=voluid;
188 fFreeParam[nmod][0]=itx;
189 fFreeParam[nmod][1]=ity;
190 fFreeParam[nmod][2]=itz;
191 fFreeParam[nmod][3]=iph;
192 fFreeParam[nmod][4]=ith;
193 fFreeParam[nmod][5]=ips;
194 nmod++;
195 }
196
197 if (strstr(st,"MODULE_VOLUID")) {
198 // to be implemented
199 }
200
201 } // end while
202
203 fNModules = nmod;
204 fNGlobal = fNModules*fgNParCh;
205
206 fclose(pfc);
207 return 0;
208}
209
210Int_t AliITSAlignMille::GetModuleIndex(const Char_t *symname) {
211 /// index from symname
212 if (!symname) return -1;
213 for (Int_t i=0; i<2198; i++) {
214 if (!strcmp(symname,AliITSgeomTGeo::GetSymName(i))) return i;
215 }
216 return -1;
217}
218
219Int_t AliITSAlignMille::GetModuleIndex(UShort_t voluid) {
220 /// index from volume ID
221 AliGeomManager::ELayerID lay = AliGeomManager::VolUIDToLayer(voluid);
222 if (lay<1|| lay>6) return -1;
223 Int_t idx=Int_t(voluid)-2048*lay;
224 if (idx>=AliGeomManager::LayerSize(lay)) return -1;
225 for (Int_t ilay=1; ilay<lay; ilay++)
226 idx += AliGeomManager::LayerSize(ilay);
227 return idx;
228}
229
230UShort_t AliITSAlignMille::GetModuleVolumeID(const Char_t *symname) {
231 /// volume ID from symname
232 if (!symname) return 0;
233
234 for (UShort_t voluid=2000; voluid<13300; voluid++) {
235 Int_t modId;
236 AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(voluid,modId);
237 if (layerId>0 && layerId<7 && modId>=0 && modId<AliGeomManager::LayerSize(layerId)) {
238 if (!strcmp(symname,AliGeomManager::SymName(layerId,modId))) return voluid;
239 }
240 }
241
242 return 0;
243}
244
245UShort_t AliITSAlignMille::GetModuleVolumeID(Int_t index) {
246 /// volume ID from index
247 if (index<0 || index>2197) return 0;
248 return GetModuleVolumeID(AliITSgeomTGeo::GetSymName(index));
249}
250
251void AliITSAlignMille::InitGeometry() {
252 /// initialize geometry
253 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
254 fGeoManager = AliGeomManager::GetGeometry();
255 if (!fGeoManager) {
256 AliInfo("Couldn't initialize geometry");
257 return;
258 }
259 // temporary align object, just use the rotation...
260 fTempAlignObj=new AliAlignObjParams(AliITSgeomTGeo::GetSymName(7),2055,0,0,0,0,0,0,kFALSE);
261}
262
263void AliITSAlignMille::Init(Int_t nGlobal, /* number of global paramers */
264 Int_t nLocal, /* number of local parameters */
265 Int_t nStdDev /* std dev cut */ )
266{
267 /// Initialization of AliMillepede. Fix parameters, define constraints ...
268 fMillepede->InitMille(nGlobal,nLocal,nStdDev,fResCut,fResCutInitial);
269 fIsMilleInit = kTRUE;
270
271 /// Fix non free parameters
272 for (Int_t i=0; i<fNModules; i++) {
273 for (Int_t j=0; j<ITSMILLE_NPARCH; j++) {
274 if (!fFreeParam[i][j]) FixParameter(i*ITSMILLE_NPARCH+j,0.0);
275 else {
276 // pepopepo: da sistemare il settaggio delle sigma...
277 Double_t parsig=0;
278 if (j<3) parsig=fParSigTranslations; // translations (0.0100 cm)
279 else parsig=fParSigRotations; // rotations (1/10 deg)
280 FixParameter(i*ITSMILLE_NPARCH+j,parsig);
281 }
282 }
283 }
284
285
286// // Set iterations
287 if (fStartFac>1) fMillepede->SetIterations(fStartFac);
288}
289
290
291void AliITSAlignMille::AddConstraint(Double_t *par, Double_t value) {
292 /// Constrain equation defined by par to value
293 if (!fIsMilleInit) {
294 AliInfo("Millepede has not been initialized!");
295 return;
296 }
297 fMillepede->SetGlobalConstraint(par, value);
298 AliInfo("Adding constraint");
299}
300
301void AliITSAlignMille::InitGlobalParameters(Double_t *par) {
302 /// Initialize global parameters with par array
303 if (!fIsMilleInit) {
304 AliInfo("Millepede has not been initialized!");
305 return;
306 }
307 fMillepede->SetGlobalParameters(par);
308 AliInfo("Init Global Parameters");
309}
310
311void AliITSAlignMille::FixParameter(Int_t iPar, Double_t value) {
312 /// Parameter iPar is encourage to vary in [-value;value].
313 /// If value == 0, parameter is fixed
314 if (!fIsMilleInit) {
315 AliInfo("Millepede has not been initialized!");
316 return;
317 }
318 fMillepede->SetParSigma(iPar, value);
319 if (value==0) AliInfo(Form("Parameter %i Fixed", iPar));
320}
321
322void AliITSAlignMille::ResetLocalEquation()
323{
324 /// Reset the derivative vectors
325 for(int i=0; i<fNLocal; i++) {
326 fLocalDerivatives[i] = 0.0;
327 }
328 for(int i=0; i<fNGlobal; i++) {
329 fGlobalDerivatives[i] = 0.0;
330 }
331}
332
333Int_t AliITSAlignMille::InitModuleParams() {
334 /// initialize geometry parameters for a given detector
335 /// for current cluster (fCluster)
336 /// fGlobalInitParam[] is set as:
337 /// [tx,ty,tz,psi,theta,phi]
338 /// (old was [tx,ty,tz,theta,psi,phi] ROOT's angles...)
339 /// *** At the moment: using Raffalele's angles definition ***
340 ///
341 /// Num of Dets: ITSMILLE_NDETELEM = fgNDetElem
342 /// Num of Par : ITSMILLE_NPARCH = fgNParCh
343 /// return 0 if success
344
345 if (!fGeoManager) {
346 AliInfo("ITS geometry not initialized!");
347 return -1;
348 }
349
350 // set the internal index (index in module list)
351 UShort_t voluid=fCluster.GetVolumeID();
352 Int_t k=fNModules-1;
353 while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
354 if (k<0) return -3;
355 fCurrentModuleInternalIndex=k;
356
357 // set the index
358 Int_t index = GetModuleIndex(AliGeomManager::SymName(voluid));
359 if (index<0) return -2;
360 fCurrentModuleIndex = index;
361
362 fModuleInitParam[0] = 0.0;
363 fModuleInitParam[1] = 0.0;
364 fModuleInitParam[2] = 0.0;
365 fModuleInitParam[3] = 0.0; // psi (X)
366 fModuleInitParam[4] = 0.0; // theta (Y)
367 fModuleInitParam[5] = 0.0; // phi (Z)
368
369 /// get global (corrected) matrix
370 // if (!AliITSgeomTGeo::GetOrigMatrix(index,*fCurrentModuleHMatrix)) return -3;
371 Double_t rott[9];
372 if (!AliITSgeomTGeo::GetRotation(index,rott)) return -3;
373 fCurrentModuleHMatrix->SetRotation(rott);
374 Double_t oLoc[3]={0,0,0};
375 if (!AliITSgeomTGeo::LocalToGlobal(index,oLoc,fCurrentModuleTranslation)) return -4;
376 fCurrentModuleHMatrix->SetTranslation(fCurrentModuleTranslation);
377
378 /// get back local coordinates
379 fMeasGlo[0] = fCluster.GetX();
380 fMeasGlo[1] = fCluster.GetY();
381 fMeasGlo[2] = fCluster.GetZ();
382 fCurrentModuleHMatrix->MasterToLocal(fMeasGlo,fMeasLoc);
383
384 // set stdev from cluster
385 TGeoHMatrix hcov;
386 Double_t hcovel[9];
387 hcovel[0]=double(fCluster.GetCov()[0]);
388 hcovel[1]=double(fCluster.GetCov()[1]);
389 hcovel[2]=double(fCluster.GetCov()[3]);
390 hcovel[3]=double(fCluster.GetCov()[1]);
391 hcovel[4]=double(fCluster.GetCov()[2]);
392 hcovel[5]=double(fCluster.GetCov()[4]);
393 hcovel[6]=double(fCluster.GetCov()[3]);
394 hcovel[7]=double(fCluster.GetCov()[4]);
395 hcovel[8]=double(fCluster.GetCov()[5]);
396 hcov.SetRotation(hcovel);
397 // now rotate in local system
398 hcov.MultiplyLeft(&fCurrentModuleHMatrix->Inverse());
399 hcov.Multiply(fCurrentModuleHMatrix);
400
401 // per i ruotati c'e' delle sigmaY che compaiono... prob
402 // e' un problema di troncamento
403 fSigmaLoc[0] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[0]));
404 fSigmaLoc[1] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[4]));
405 fSigmaLoc[2] = TMath::Sqrt(TMath::Abs(hcov.GetRotationMatrix()[8]));
406
407 AliDebug(2,Form("Setting StDev from CovMat : fSigmaLocX=%f fSigmaLocY=%f fSigmaLocZ=%f \n",fSigmaLoc[0] ,fSigmaLoc[1] ,fSigmaLoc[2] ));
408
409 return 0;
410}
411
412void AliITSAlignMille::SetCurrentModule(Int_t index) {
413 ///
414 UShort_t voluid=GetModuleVolumeID(index);
415 if (voluid) {
416 fCluster.SetVolumeID(voluid);
417 fCluster.SetXYZ(0,0,0);
418 InitModuleParams();
419 }
420}
421
422void AliITSAlignMille::Print() {
423 ///
424 printf("*** AliMillepede for ITS ***\n");
425 printf(" number of defined modules: %d\n",fNModules);
426 if (fGeoManager)
427 printf(" geometry loaded from %s\n",fGeometryFileName.Data());
428 else
429 printf(" geometry not loaded\n");
430 if (fUseLocalShifts)
431 printf(" Alignment shifts will be computed in LOCAL RS\n");
432 else
433 printf(" Alignment shifts will be computed in GLOBAL RS\n");
434 printf(" Millepede configuration parameters:\n");
435 printf(" init parsig for translations : %.4f\n",fParSigTranslations);
436 printf(" init parsig for rotations : %.4f\n",fParSigRotations);
437 printf(" init value for chi2 cut : %.4f\n",fStartFac);
438 printf(" first iteration cut value : %.4f\n",fResCutInitial);
439 printf(" other iterations cut value : %.4f\n",fResCut);
440 printf(" number of stddev for chi2 cut : %d\n",fNStdDev);
441
442 printf("List of defined modules:\n");
443 printf(" intidx\tindex\tvoluid\tsymname\n");
444 for (int i=0; i<fNModules; i++)
445 printf(" %d\t%d\t%d\t%s\n",i,fModuleIndex[i],fModuleVolumeID[i],AliITSgeomTGeo::GetSymName(fModuleIndex[i]));
446}
447
448void AliITSAlignMille::PrintCurrentModuleInfo() {
449 ///
450 if (fCurrentModuleIndex<0 || fCurrentModuleIndex>2197) return;
451 UShort_t voluid=fModuleVolumeID[fCurrentModuleInternalIndex];
452 printf("Current module: index=%d voluid=%d\n",fCurrentModuleIndex,voluid);
453 printf(" symname:%s\n",AliGeomManager::SymName(voluid));
454 printf(" TGeoHMatrix: \n");
455 fCurrentModuleHMatrix->Print();
456}
457
458
459void AliITSAlignMille::InitTrackParams(int meth) {
460 /// initialize local parameters with different methods
461 /// for current track (fTrack)
462
463 switch (meth) {
464 case 1: // simple linear interpolation
465 // get local starting parameters (to be substituted by ESD track parms)
466 // local parms (fLocalInitParam[]) are:
467 // [0] = global x coord. of straight line intersection at y=0 plane
468 // [1] = global z coord. of straight line intersection at y=0 plane
469 // [2] = px/py
470 // [3] = pz/py
471
472 // test #1: linear fit in x(y) and z(y)
473 Int_t npts = fTrack->GetNPoints();
474
475 TF1 *f1=new TF1("f1","[0]+x*[1]",-50,50);
476
477 TGraph *g=new TGraph(npts,fTrack->GetY(),fTrack->GetX());
478 g->Fit(f1,"RNQ");
479 fLocalInitParam[0] = f1->GetParameter(0);
480 fLocalInitParam[2] = f1->GetParameter(1);
481 AliDebug(2,Form("X = p0gx + ugx*Y : p0gx = %f +- %f ugx = %f +- %f\n",fLocalInitParam[0],f1->GetParError(0),fLocalInitParam[2],f1->GetParError(1)));
482 delete g; g=NULL;
483
484 g=new TGraph(npts,fTrack->GetY(),fTrack->GetZ());
485 g->Fit(f1,"RNQ");
486 fLocalInitParam[1] = f1->GetParameter(0);
487 fLocalInitParam[3] = f1->GetParameter(1);
488 AliDebug(2,Form("Z = p0gz + ugz*Y : p0gz=%f ugz=%f\n",fLocalInitParam[1],fLocalInitParam[3]));
489 delete g;
490 delete f1;
491
492 break;
493 }
494
495}
496Bool_t AliITSAlignMille::CheckVolumeID(UShort_t voluid) const
497{
498 ///
499 Int_t k=fNModules-1;
500 while (k>=0 && !(voluid==fModuleVolumeID[k]) ) k--;
501 //printf("selected element with voluid=%d : %d\n",voluid,k);
502 if (k>=0) return kTRUE;
503 return kFALSE;
504}
505
506Int_t AliITSAlignMille::CheckCurrentTrack() {
507 /// checks if AliTrackPoints belongs to defined modules
508 /// return number of good poins
509 /// return 0 if not enough points
510
511 Int_t npts = fTrack->GetNPoints();
512 Int_t ngoodpts=0;
513 // debug points
514 for (int j=0; j<npts; j++) {
515 UShort_t voluid = fTrack->GetVolumeID()[j];
516 if (CheckVolumeID(voluid)) {
517 ngoodpts++;
518 }
519 }
520 // pepo da controllare...
521 if (ngoodpts<fMinNPtsPerTrack) return 0;
522
523 return ngoodpts;
524}
525
526Int_t AliITSAlignMille::ProcessTrack(AliTrackPointArray *track) {
527 /// Process track; Loop over hits and set local equations
528 /// here 'track' is a AliTrackPointArray
529 /// return 0 if success;
530
531 if (!fIsMilleInit) {
532 AliInfo("Millepede has not been initialized!");
533 return -1;
534 }
535
536 Int_t npts = track->GetNPoints();
537 AliDebug(2,Form("\n*** Processing track with %d points ***\n",npts));
538 fTrack = track;
539
540 // checks if there are enough good points
541 if (!CheckCurrentTrack()) {
542 AliInfo("Track with not enough good points - will not be used...");
543 return -1;
544 }
545
546 // set local starting parameters (to be substituted by ESD track parms)
547 // local parms (fLocalInitParam[]) are:
548 // [0] = global x coord. of straight line intersection at y=0 plane
549 // [1] = global z coord. of straight line intersection at y=0 plane
550 // [2] = px/py
551 // [3] = pz/py
552 InitTrackParams(1);
553
554 for (Int_t ipt=0; ipt<npts; ipt++) {
555 fTrack->GetPoint(fCluster,ipt);
556 if (!CheckVolumeID(fCluster.GetVolumeID())) continue;
557
558 // set geometry parameters for the the current module
559 AliDebug(2,Form("\n--- processing point %d --- \n",ipt));
560 if (InitModuleParams()) continue;
561
562 AliDebug(2,Form(" VolID=%d Index=%d InternalIdx=%d symname=%s\n", track->GetVolumeID()[ipt], fCurrentModuleIndex ,fCurrentModuleInternalIndex, AliGeomManager::SymName(track->GetVolumeID()[ipt]) ));
563 AliDebug(2,Form(" Point = ( %f , %f , %f ) \n",track->GetX()[ipt],track->GetY()[ipt],track->GetZ()[ipt]));
564
565 if (SetLocalEquations()) return -1;
566
567 } // end loop over points
568
569 return 0;
570}
571
572Int_t AliITSAlignMille::CalcIntersectionPoint(Double_t *lpar, Double_t *gpar) {
573 /// calculate track intersection point in local coordinates
574 /// according with a given set of parameters (local(4) and global(6))
575 /// and fill fPintLoc/Glo
576 /// local are: pgx0, pgz0, ugx0, ugz0
577 /// global are: tx,ty,tz,psi,theta,phi (Raff's delta angles in deg.)
578 /// return 0 if success
579
580 AliDebug(3,Form("lpar = %g %g %g %g \ngpar= %g %g %g %g %g %g\n",lpar[0],lpar[1],lpar[2],lpar[3],gpar[0],gpar[1],gpar[2],gpar[3],gpar[4],gpar[5]));
581
582 // vector of initial straight line direction in glob. coord
583 // ATTENZIONE: forse va rinormalizzato tutto...
584 Double_t v0g[3];
585 //Double_t
586 v0g[0]=lpar[2];
587 v0g[1]=1.0;
588 v0g[2]=lpar[3];
589
590 // intercept in yg=0 plane in glob coord
591 Double_t p0g[3];
592 p0g[0]=lpar[0];
593 p0g[1]=0.0;
594 p0g[2]=lpar[1];
595
596
597 // prepare the TGeoHMatrix
598 Double_t tr[3],ang[3];
599 //Double_t rad2deg=180./TMath::Pi();
600 if (fUseLocalShifts) { // just Delta matrix
601 tr[0]=gpar[0];
602 tr[1]=gpar[1];
603 tr[2]=gpar[2];
604 ang[0]=gpar[3]; // psi (X)
605 ang[1]=gpar[4]; // theta (Y)
606 ang[2]=gpar[5]; // phi (Z)
607 }
608 else { // total matrix with shifted parameter
609 AliInfo("global shifts not implemented yet!");
610 return -1;
611 }
612
613 //printf("fTempRot = 0x%x - ang = %g %g %g \n",fTempRot,gpar[5]*rad2deg,gpar[3]*rad2deg,gpar[4]*rad2deg);
614
615 fTempAlignObj->SetRotation(ang[0],ang[1],ang[2]);
616 AliDebug(3,Form("Delta angles: psi=%f theta=%f phi=%f",ang[0],ang[1],ang[2]));
617 TGeoHMatrix hm;
618 fTempAlignObj->GetMatrix(hm);
619 fTempHMat->SetRotation(hm.GetRotationMatrix());
620 fTempHMat->SetTranslation(tr);
621
622 // in this case the gpar[] array contains only shifts
623 // and fInitModuleParam[] are set to 0
624 // fCurrentModuleHMatrix is then modified as fCurrentHM*fTempHM
625 if (fUseLocalShifts)
626 fTempHMat->MultiplyLeft(fCurrentModuleHMatrix);
627
628 // same in local coord.
629 Double_t p0l[3],v0l[3];
630 fTempHMat->MasterToLocalVect(v0g,v0l);
631 fTempHMat->MasterToLocal(p0g,p0l);
632
633 if (TMath::Abs(v0l[1])<1e-15) {
634 AliInfo("Track Y direction in local frame is zero! Cannot proceed...");
635 return -1;
636 }
637
638 // local intersection point
639 fPintLoc[0] = p0l[0] - (v0l[0]/v0l[1])*p0l[1];
640 fPintLoc[1] = 0;
641 fPintLoc[2] = p0l[2] - (v0l[2]/v0l[1])*p0l[1];
642
643 // global intersection point
644 fTempHMat->LocalToMaster(fPintLoc,fPintGlo);
645 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]));
646
647 return 0;
648}
649
650Int_t AliITSAlignMille::CalcDerivatives(Int_t paridx, Bool_t islpar) {
651 /// calculate numerically (ROOT's style) the derivatives for
652 /// local X intersection and local Z intersection
653 /// parlist: local (islpar=kTRUE) pgx0, pgz0, ugx0, ugz0
654 /// global (islpar=kFALSE) tx, ty, tz, psi, theta, phi (Raf's angles in deg)
655 /// return 0 if success
656
657 // copy initial parameters
658 Double_t lpar[ITSMILLE_NLOCAL];
659 Double_t gpar[ITSMILLE_NPARCH];
660 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) lpar[i]=fLocalInitParam[i];
661 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) gpar[i]=fModuleInitParam[i];
662
663 // pepopepo
664 // trial with fixed dpar...
665 Double_t dpar=0.0;
666 if (islpar) {
667 //dpar=fLocalInitParam[paridx]*0.001;
668 // set minimum dpar
669 if (paridx<2) dpar=1.0e-4; // translations
670 else dpar=1.0e-6; // direction
671 }
672 else {
673 //dpar=fModuleInitParam[paridx]*0.001;
674 if (paridx<3) dpar=1.0e-4; // translations
675 else dpar=1.0e-2; // angles
676 }
677 AliDebug(3,Form("\n+++ automatic dpar=%g\n",dpar));
678 if (fDeltaPar) dpar=fDeltaPar;
679 AliDebug(3,Form("+++ using dpar=%g\n\n",dpar));
680
681 // calculate derivative ROOT's like:
682 // using f(x+h),f(x-h),f(x+h/2),f(x-h2)...
683 Double_t pintl1[3]; // f(x-h)
684 Double_t pintl2[3]; // f(x-h/2)
685 Double_t pintl3[3]; // f(x+h/2)
686 Double_t pintl4[3]; // f(x+h)
687
688 // first values
689 if (islpar) lpar[paridx] -= dpar;
690 else gpar[paridx] -= dpar;
691 if (CalcIntersectionPoint(lpar, gpar)) return -2;
692 for (Int_t i=0; i<3; i++) pintl1[i]=fPintLoc[i];
693
694 // second values
695 if (islpar) lpar[paridx] += dpar/2;
696 else gpar[paridx] += dpar/2;
697 if (CalcIntersectionPoint(lpar, gpar)) return -2;
698 for (Int_t i=0; i<3; i++) pintl2[i]=fPintLoc[i];
699
700 // third values
701 if (islpar) lpar[paridx] += dpar;
702 else gpar[paridx] += dpar;
703 if (CalcIntersectionPoint(lpar, gpar)) return -2;
704 for (Int_t i=0; i<3; i++) pintl3[i]=fPintLoc[i];
705
706 // fourth values
707 if (islpar) lpar[paridx] += dpar/2;
708 else gpar[paridx] += dpar/2;
709 if (CalcIntersectionPoint(lpar, gpar)) return -2;
710 for (Int_t i=0; i<3; i++) pintl4[i]=fPintLoc[i];
711
712 Double_t h2 = 1./(2.*dpar);
713 Double_t d0 = pintl4[0]-pintl1[0];
714 Double_t d2 = 2.*(pintl3[0]-pintl2[0]);
715 fDerivativeXLoc = h2*(4*d2 - d0)/3.;
716 if (TMath::Abs(fDerivativeXLoc) < 1.0e-9) fDerivativeXLoc=0.0;
717
718 d0 = pintl4[2]-pintl1[2];
719 d2 = 2.*(pintl3[2]-pintl2[2]);
720 fDerivativeZLoc = h2*(4*d2 - d0)/3.;
721 if (TMath::Abs(fDerivativeZLoc) < 1.0e-9) fDerivativeZLoc=0.0;
722
723 AliDebug(3,Form("\n+++ derivatives +++ \n"));
724 AliDebug(3,Form("+++ dXLoc/dpar = %g +++\n",fDerivativeXLoc));
725 AliDebug(3,Form("+++ dZLoc/dpar = %g +++\n\n",fDerivativeZLoc));
726
727 return 0;
728}
729
730Int_t AliITSAlignMille::SetLocalEquations() {
731 /// Define local equation for current track and cluster in x coor.
732 /// return 0 if success
733
734 // store first interaction point
735 CalcIntersectionPoint(fLocalInitParam, fModuleInitParam);
736 for (Int_t i=0; i<3; i++) fPintLoc0[i]=fPintLoc[i];
737
738 // calculate local derivatives numerically
739 Double_t dXdL[ITSMILLE_NLOCAL],dZdL[ITSMILLE_NLOCAL];
740 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++) {
741 if (CalcDerivatives(i,kTRUE)) return -1;
742 dXdL[i]=fDerivativeXLoc;
743 dZdL[i]=fDerivativeZLoc;
744 }
745
746 Double_t dXdG[ITSMILLE_NPARCH],dZdG[ITSMILLE_NPARCH];
747 for (Int_t i=0; i<ITSMILLE_NPARCH; i++) {
748 if (CalcDerivatives(i,kFALSE)) return -1;
749 dXdG[i]=fDerivativeXLoc;
750 dZdG[i]=fDerivativeZLoc;
751 }
752
753 AliDebug(2,Form("\n***************\n"));
754 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
755 AliDebug(2,Form("Local parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdL[i],dZdL[i]));
756 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
757 AliDebug(2,Form("Global parameter %d - dXdpar = %g - dZdpar = %g\n",i,dXdG[i],dZdG[i]));
758 AliDebug(2,Form("\n***************\n"));
759
760
761 AliDebug(2,Form("setting local equation X with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]));
762 // set equation for Xloc coordinate
763 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
764 SetLocalDerivative(i,dXdL[i]);
765 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
766 SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dXdG[i]);
767 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[0]-fPintLoc0[0]), fSigmaLoc[0]);
768
769
770 AliDebug(2,Form("setting local equation Z with fMeas=%.6f and fSigma=%.6f",(fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]));
771 // set equation for Zloc coordinate
772 for (Int_t i=0; i<ITSMILLE_NLOCAL; i++)
773 SetLocalDerivative(i,dZdL[i]);
774 for (Int_t i=0; i<ITSMILLE_NPARCH; i++)
775 SetGlobalDerivative(fCurrentModuleInternalIndex*ITSMILLE_NPARCH+i,dZdG[i]);
776 fMillepede->SetLocalEquation(fGlobalDerivatives, fLocalDerivatives, (fMeasLoc[2]-fPintLoc0[2]), fSigmaLoc[2]);
777
778 return 0;
779}
780
781
782void AliITSAlignMille::LocalFit(Int_t iTrack, Double_t *lTrackParam, Int_t lSingleFit) {
783 /// Call local fit for this track
784 if (!fIsMilleInit) {
785 AliInfo("Millepede has not been initialized!");
786 return;
787 }
788 Int_t iRes = fMillepede->LocalFit(iTrack,lTrackParam,lSingleFit);
789 AliDebug(2,Form("iRes = %d",iRes));
790 if (iRes && !lSingleFit) {
791 fMillepede->SetNLocalEquations(fMillepede->GetNLocalEquations()+1);
792 }
793}
794
795void AliITSAlignMille::GlobalFit(Double_t *parameters,Double_t *errors,Double_t *pulls) {
796 /// Call global fit; Global parameters are stored in parameters
797 if (!fIsMilleInit) {
798 AliInfo("Millepede has not been initialized!");
799 return;
800 }
801 fMillepede->GlobalFit(parameters,errors,pulls);
802 AliInfo("Done fitting global parameters!");
803}
804
805Double_t AliITSAlignMille::GetParError(Int_t iPar) {
806 /// Get error of parameter iPar
807 if (!fIsMilleInit) {
808 AliInfo("Millepede has not been initialized!");
809 return 0;
810 }
811 Double_t lErr = fMillepede->GetParError(iPar);
812 return lErr;
813}
814
815void AliITSAlignMille::PrintGlobalParameters() {
816 /// Print global parameters
817 if (!fIsMilleInit) {
818 AliInfo("Millepede has not been initialized!");
819 return;
820 }
821 fMillepede->PrintGlobalParameters();
822}
823
824// //_________________________________________________________________________
825