]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliAlignmentTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Change from Int_t/SHort_t to smaller data types, the change from Float_t/Double_t...
[u/mrichter/AliRoot.git] / STEER / AliAlignmentTracks.cxx
CommitLineData
98937d93 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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//-----------------------------------------------------------------
17// Implementation of the alignment steering class
18// It provides an access to the track space points
19// written along the esd tracks. The class enables
20// the user to plug any track fitter (deriving from
21// AliTrackFitter class) and minimization fo the
22// track residual sums (deriving from the AliTrackResiduals).
23//-----------------------------------------------------------------
24
25#include <TChain.h>
bfbd5665 26#include <TFile.h>
d6116b3f 27#include <TVector3.h>
98937d93 28
29#include "AliAlignmentTracks.h"
30#include "AliTrackPointArray.h"
90dbf5fb 31#include "AliAlignObjParams.h"
98937d93 32#include "AliTrackFitterRieman.h"
33#include "AliTrackResidualsChi2.h"
9a1304f0 34#include "AliESDEvent.h"
98937d93 35#include "AliLog.h"
36
37ClassImp(AliAlignmentTracks)
38
39//______________________________________________________________________________
40AliAlignmentTracks::AliAlignmentTracks():
41 fESDChain(0),
42 fPointsFilename("AliTrackPoints.root"),
43 fPointsFile(0),
44 fPointsTree(0),
45 fLastIndex(0),
46 fArrayIndex(0),
47 fIsIndexBuilt(kFALSE),
75e3794b 48 fAlignObjs(0),
46ae650f 49 fMisalignObjs(0),
98937d93 50 fTrackFitter(0),
91264393 51 fMinimizer(0),
52 fDoUpdate(kTRUE)
98937d93 53{
54 // Default constructor
55 InitIndex();
56 InitAlignObjs();
57}
58
59//______________________________________________________________________________
60AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
61 fESDChain(esdchain),
62 fPointsFilename("AliTrackPoints.root"),
63 fPointsFile(0),
64 fPointsTree(0),
65 fLastIndex(0),
66 fArrayIndex(0),
67 fIsIndexBuilt(kFALSE),
75e3794b 68 fAlignObjs(0),
46ae650f 69 fMisalignObjs(0),
98937d93 70 fTrackFitter(0),
91264393 71 fMinimizer(0),
72 fDoUpdate(kTRUE)
98937d93 73{
74 // Constructor in the case
75 // the user provides an already
76 // built TChain with ESD trees
77 InitIndex();
78 InitAlignObjs();
79}
80
81
82//______________________________________________________________________________
83AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
75e3794b 84 fESDChain(new TChain(esdtreename)),
98937d93 85 fPointsFilename("AliTrackPoints.root"),
86 fPointsFile(0),
87 fPointsTree(0),
88 fLastIndex(0),
89 fArrayIndex(0),
90 fIsIndexBuilt(kFALSE),
75e3794b 91 fAlignObjs(0),
46ae650f 92 fMisalignObjs(0),
98937d93 93 fTrackFitter(0),
91264393 94 fMinimizer(0),
95 fDoUpdate(kTRUE)
98937d93 96{
97 // Constructor in the case
98 // the user provides a single ESD file
99 // or a directory containing ESD files
98937d93 100 fESDChain->Add(esdfilename);
101
102 InitIndex();
103 InitAlignObjs();
104}
105
98937d93 106
107//______________________________________________________________________________
108AliAlignmentTracks::~AliAlignmentTracks()
109{
110 // Destructor
111 if (fESDChain) delete fESDChain;
112
113 DeleteIndex();
114 DeleteAlignObjs();
115
116 delete fTrackFitter;
117 delete fMinimizer;
118
119 if (fPointsFile) fPointsFile->Close();
120}
121
122//______________________________________________________________________________
123void AliAlignmentTracks::AddESD(TChain *esdchain)
124{
125 // Add a chain with ESD files
126 if (fESDChain)
127 fESDChain->Add(esdchain);
128 else
129 fESDChain = esdchain;
130}
131
132//______________________________________________________________________________
133void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
134{
135 // Add a single file or
136 // a directory to the chain
137 // with the ESD files
138 if (fESDChain)
139 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
140 else {
141 fESDChain = new TChain(esdtreename);
142 fESDChain->Add(esdfilename);
143 }
144}
145
babbb915 146//________________________________________________________________________
babbb915 147void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
148 Int_t minITSpts,
149 Bool_t cuts,
d6116b3f 150 Float_t minAngleWrtITSModulePlanes,
babbb915 151 Float_t minMom,Float_t maxMom,
152 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
153 Float_t minSinTheta,Float_t maxSinTheta)
98937d93 154{
155 // Analyzes and filters ESD tracks
156 // Stores the selected track space points
157 // into the output file
158
159 if (!fESDChain) return;
160
9a1304f0 161 AliESDEvent *esd = new AliESDEvent();
162 esd->ReadFromTree(fESDChain);
cf0f66c2 163 AliESDfriend *esdf = 0;
164 fESDChain->SetBranchStatus("ESDfriend*",1);
165 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
98937d93 166
167 // Open the output file
168 if (fPointsFilename.Data() == "") {
169 AliWarning("Incorrect output filename!");
170 return;
171 }
172
173 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
174 if (!pointsFile || !pointsFile->IsOpen()) {
175 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
176 return;
177 }
178
179 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
15e85efa 180 const AliTrackPointArray *array = 0;
babbb915 181 AliTrackPointArray *array2 = 0;
182 if(onlyITS) { // only ITS AliTrackPoints
183 pointsTree->Branch("SP","AliTrackPointArray", &array2);
184 } else {
185 pointsTree->Branch("SP","AliTrackPointArray", &array);
186 }
187
98937d93 188 Int_t ievent = 0;
189 while (fESDChain->GetEntry(ievent++)) {
190 if (!esd) break;
cf0f66c2 191
192 esd->SetESDfriend(esdf); //Attach the friend to the ESD
193
98937d93 194 Int_t ntracks = esd->GetNumberOfTracks();
195 for (Int_t itrack=0; itrack < ntracks; itrack++) {
196 AliESDtrack * track = esd->GetTrack(itrack);
197 if (!track) continue;
98937d93 198
babbb915 199 if(track->GetNcls(0) < minITSpts) continue;
200 if(cuts) {
201 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
202 Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
203 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
204 Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
205 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
206 }
98937d93 207
babbb915 208 AliTrackPoint point;
98937d93 209 array = track->GetTrackPointArray();
babbb915 210
211 if(onlyITS) {
212 array2 = new AliTrackPointArray(track->GetNcls(0));
d6116b3f 213 Bool_t smallAngleWrtITSModulePlanes=kFALSE;
babbb915 214 Int_t ipt,volId,modId,layerId;
215 Int_t jpt=0;
216 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
217 array->GetPoint(point,ipt);
218 volId = point.GetVolumeID();
219 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
220 if(layerId<7) {
221 array2->AddPoint(jpt,&point);
222 jpt++;
223 }
d6116b3f 224 // check minAngleWrtITSModulePlanes
225 if(cuts) {
226 Double_t p[3]; track->GetDirection(p);
227 TVector3 pvec(p[0],p[1],p[2]);
228 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
229 TVector3 normvec(rot[1],-rot[0],0.);
230 Double_t angle = pvec.Angle(normvec);
231 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
232 angle = 0.5*TMath::Pi()-angle;
233 if(angle<minAngleWrtITSModulePlanes) smallAngleWrtITSModulePlanes=kTRUE;
234 }
babbb915 235 }
d6116b3f 236 if(smallAngleWrtITSModulePlanes) continue;
babbb915 237 } // end if(onlyITS)
238
98937d93 239 pointsTree->Fill();
240 }
241 }
242
243 if (!pointsTree->Write()) {
244 AliWarning("Can't write the tree with track point arrays!");
245 return;
246 }
247
248 pointsFile->Close();
d6116b3f 249
250 return;
98937d93 251}
252
d6116b3f 253//_____________________________________________________________________________
babbb915 254void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
9a1304f0 255 Int_t minITSpts,Float_t maxMatchingAngle,
d6116b3f 256 Bool_t cuts,
257 Float_t minAngleWrtITSModulePlanes,
babbb915 258 Float_t minMom,Float_t maxMom,
259 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
260 Float_t minSinTheta,Float_t maxSinTheta)
261{
262 // Analyzes and filters ESD tracks
263 // Merges inward and outward tracks in one single track
264 // Stores the selected track space points
265 // into the output file
266
267 if (!fESDChain) return;
268
9a1304f0 269 AliESDEvent *esd = new AliESDEvent();
270 esd->ReadFromTree(fESDChain);
babbb915 271 AliESDfriend *esdf = 0;
272 fESDChain->SetBranchStatus("ESDfriend*",1);
273 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
274
275 // Open the output file
276 if (fPointsFilename.Data() == "") {
277 AliWarning("Incorrect output filename!");
278 return;
279 }
280
281 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
282 if (!pointsFile || !pointsFile->IsOpen()) {
283 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
284 return;
285 }
286
287 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
288 const AliTrackPointArray *array = 0;
289 AliTrackPointArray *array2 = 0;
290 pointsTree->Branch("SP","AliTrackPointArray", &array2);
291
292 Int_t ievent = 0;
293 while (fESDChain->GetEntry(ievent++)) {
294 if (!esd) break;
295
296 esd->SetESDfriend(esdf); //Attach the friend to the ESD
297
298 Int_t ntracks = esd->GetNumberOfTracks();
9a1304f0 299 if(ntracks<2) continue;
300 Int_t *goodtracksArray = new Int_t[ntracks];
301 Float_t *phiArray = new Float_t[ntracks];
302 Float_t *thetaArray = new Float_t[ntracks];
babbb915 303 Int_t ngt=0;
babbb915 304 for (Int_t itrack=0; itrack < ntracks; itrack++) {
305 AliESDtrack * track = esd->GetTrack(itrack);
306 if (!track) continue;
307
308 if(track->GetNcls(0) < minITSpts) continue;
9a1304f0 309 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
310 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
babbb915 311 if(cuts) {
312 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
9a1304f0 313 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
babbb915 314 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
9a1304f0 315 Float_t sintheta = TMath::Sin(theta);
babbb915 316 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
317 }
9a1304f0 318 goodtracksArray[ngt]=itrack;
319 phiArray[ngt]=phi;
320 thetaArray[ngt]=theta;
babbb915 321 ngt++;
322 }
323
9a1304f0 324 if(ngt<2) {
325 delete [] goodtracksArray; goodtracksArray=0;
326 delete [] phiArray; phiArray=0;
327 delete [] thetaArray; thetaArray=0;
328 continue;
329 }
330
331 // check matching of the two tracks from the muon
332 Float_t min = 10000000.;
333 Int_t good1 = -1, good2 = -1;
334 for(Int_t itr1=0; itr1<ngt-1; itr1++) {
335 for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
336 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
337 if(deltatheta>maxMatchingAngle) continue;
338 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
339 if(deltaphi>maxMatchingAngle) continue;
d6116b3f 340 //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
9a1304f0 341 if(deltatheta+deltaphi<min) {
342 min=deltatheta+deltaphi;
343 good1 = goodtracksArray[itr1];
344 good2 = goodtracksArray[itr2];
345 }
346 }
347 }
348
349 delete [] goodtracksArray; goodtracksArray=0;
350 delete [] phiArray; phiArray=0;
351 delete [] thetaArray; thetaArray=0;
352
353 if(good1<0) continue;
babbb915 354
9a1304f0 355 AliESDtrack * track1 = esd->GetTrack(good1);
356 AliESDtrack * track2 = esd->GetTrack(good2);
babbb915 357
358 Int_t ntotpts;
359 if(onlyITS) {
360 ntotpts = track1->GetNcls(0)+track2->GetNcls(0);
361 } else {
362 ntotpts = track1->GetTrackPointArray()->GetNPoints()+
363 track2->GetTrackPointArray()->GetNPoints();
364 }
365 array2 = new AliTrackPointArray(ntotpts);
366 AliTrackPoint point;
367 Int_t ipt,volId,modId,layerId;
368 Int_t jpt=0;
d6116b3f 369 Bool_t smallAngleWrtITSModulePlanes=kFALSE;
babbb915 370 array = track1->GetTrackPointArray();
371 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
372 array->GetPoint(point,ipt);
373 if(onlyITS) {
374 volId = point.GetVolumeID();
375 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
376 if(layerId>6) continue;
d6116b3f 377 // check minAngleWrtITSModulePlanes
378 if(cuts) {
379 Double_t p[3]; track1->GetDirection(p);
380 TVector3 pvec(p[0],p[1],p[2]);
381 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
382 TVector3 normvec(rot[1],-rot[0],0.);
383 Double_t angle = pvec.Angle(normvec);
384 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
385 angle = 0.5*TMath::Pi()-angle;
386 if(angle<minAngleWrtITSModulePlanes) smallAngleWrtITSModulePlanes=kTRUE;
387 }
babbb915 388 }
389 array2->AddPoint(jpt,&point);
390 jpt++;
391 }
392 array = track2->GetTrackPointArray();
393 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
394 array->GetPoint(point,ipt);
395 if(onlyITS) {
396 volId = point.GetVolumeID();
397 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
398 if(layerId>6) continue;
d6116b3f 399 // check minAngleWrtITSModulePlanes
400 if(cuts) {
401 Double_t p[3]; track2->GetDirection(p);
402 TVector3 pvec(p[0],p[1],p[2]);
403 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
404 TVector3 normvec(rot[1],-rot[0],0.);
405 Double_t angle = pvec.Angle(normvec);
406 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
407 angle = 0.5*TMath::Pi()-angle;
408 if(angle<minAngleWrtITSModulePlanes) smallAngleWrtITSModulePlanes=kTRUE;
409 }
babbb915 410 }
411 array2->AddPoint(jpt,&point);
412 jpt++;
413 }
414
d6116b3f 415 if(smallAngleWrtITSModulePlanes) continue;
babbb915 416 pointsTree->Fill();
417 }
418
419 if (!pointsTree->Write()) {
420 AliWarning("Can't write the tree with track point arrays!");
421 return;
422 }
423
424 pointsFile->Close();
d6116b3f 425 return;
babbb915 426}
427
98937d93 428//______________________________________________________________________________
429void AliAlignmentTracks::ProcessESD(TSelector *selector)
430{
431 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
432}
433
434//______________________________________________________________________________
435void AliAlignmentTracks::BuildIndex()
436{
437 // Build index of points tree entries
438 // Used for access based on the volume IDs
46ae650f 439 if (fIsIndexBuilt) return;
440
441 fIsIndexBuilt = kTRUE;
98937d93 442
33bebad6 443 // Dummy object is created in order
444 // to initialize the volume paths
90dbf5fb 445 AliAlignObjParams alobj;
33bebad6 446
7fa594d2 447 fPointsFile = TFile::Open(fPointsFilename);
98937d93 448 if (!fPointsFile || !fPointsFile->IsOpen()) {
449 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
450 return;
451 }
452
453 // AliTrackPointArray* array = new AliTrackPointArray;
454 AliTrackPointArray* array = 0;
46ae650f 455 fPointsTree = (TTree*) fPointsFile->Get("spTree");
456 if (!fPointsTree) {
98937d93 457 AliWarning("No pointsTree found!");
458 return;
459 }
46ae650f 460 fPointsTree->SetBranchAddress("SP", &array);
98937d93 461
7fa594d2 462 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
98937d93 463 for (Int_t iArray = 0; iArray < nArrays; iArray++)
464 {
46ae650f 465 fPointsTree->GetEvent(iArray);
98937d93 466 if (!array) continue;
467 for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
468 UShort_t volId = array->GetVolumeID()[ipoint];
33bebad6 469 // check if the volId is valid
25be1e5c 470 if (!AliGeomManager::SymName(volId)) {
b760c02e 471 AliError(Form("The volume id %d has no default volume name !",
6baf24ce 472 volId));
473 continue;
474 }
98937d93 475 Int_t modId;
25be1e5c 476 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
477 - AliGeomManager::kFirstLayer;
98937d93 478 if (!fArrayIndex[layerId][modId]) {
479 //first entry for this volume
480 fArrayIndex[layerId][modId] = new TArrayI(1000);
481 }
482 else {
483 Int_t size = fArrayIndex[layerId][modId]->GetSize();
484 // If needed allocate new size
485 if (fLastIndex[layerId][modId] >= size)
486 fArrayIndex[layerId][modId]->Set(size + 1000);
487 }
488
489 // Check if the index is already filled
490 Bool_t fillIndex = kTRUE;
491 if (fLastIndex[layerId][modId] != 0) {
492 if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
493 fillIndex = kFALSE;
494 }
495 // Fill the index array and store last filled index
496 if (fillIndex) {
497 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
498 fLastIndex[layerId][modId]++;
499 }
500 }
501 }
502
503}
504
98937d93 505//______________________________________________________________________________
506void AliAlignmentTracks::InitIndex()
507{
508 // Initialize the index arrays
25be1e5c 509 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
98937d93 510 fLastIndex = new Int_t*[nLayers];
511 fArrayIndex = new TArrayI**[nLayers];
25be1e5c 512 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
513 fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
514 fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
515 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
98937d93 516 fLastIndex[iLayer][iModule] = 0;
517 fArrayIndex[iLayer][iModule] = 0;
518 }
519 }
520}
521
522//______________________________________________________________________________
523void AliAlignmentTracks::ResetIndex()
524{
525 // Reset the value of the last filled index
526 // Do not realocate memory
46ae650f 527
528 fIsIndexBuilt = kFALSE;
529
25be1e5c 530 for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
531 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
98937d93 532 fLastIndex[iLayer][iModule] = 0;
533 }
534 }
535}
536
537//______________________________________________________________________________
538void AliAlignmentTracks::DeleteIndex()
539{
540 // Delete the index arrays
541 // Called by the destructor
25be1e5c 542 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
543 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
98937d93 544 if (fArrayIndex[iLayer][iModule]) {
545 delete fArrayIndex[iLayer][iModule];
546 fArrayIndex[iLayer][iModule] = 0;
547 }
548 }
549 delete [] fLastIndex[iLayer];
550 delete [] fArrayIndex[iLayer];
551 }
552 delete [] fLastIndex;
553 delete [] fArrayIndex;
554}
555
556//______________________________________________________________________________
46ae650f 557Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
98937d93 558{
559 // Read alignment object from a file
560 // To be replaced by a call to CDB
46ae650f 561 AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
98937d93 562
563 return kFALSE;
564}
565
566//______________________________________________________________________________
567void AliAlignmentTracks::InitAlignObjs()
568{
569 // Initialize the alignment objects array
25be1e5c 570 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
98937d93 571 fAlignObjs = new AliAlignObj**[nLayers];
25be1e5c 572 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
573 fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
574 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
575 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
90dbf5fb 576 fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
46ae650f 577 }
98937d93 578 }
579}
580
581//______________________________________________________________________________
582void AliAlignmentTracks::ResetAlignObjs()
583{
584 // Reset the alignment objects array
25be1e5c 585 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
586 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
98937d93 587 fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
588 }
589}
590
591//______________________________________________________________________________
592void AliAlignmentTracks::DeleteAlignObjs()
593{
594 // Delete the alignment objects array
25be1e5c 595 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
596 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
98937d93 597 if (fAlignObjs[iLayer][iModule])
598 delete fAlignObjs[iLayer][iModule];
599 delete [] fAlignObjs[iLayer];
600 }
601 delete [] fAlignObjs;
46ae650f 602 fAlignObjs = 0;
98937d93 603}
604
25be1e5c 605void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
606 AliGeomManager::ELayerID lastLayer,
607 AliGeomManager::ELayerID layerRangeMin,
608 AliGeomManager::ELayerID layerRangeMax,
cc345ce3 609 Int_t iterations)
98937d93 610{
cc345ce3 611 // Align detector volumes within
612 // a given layer range
613 // (could be whole detector).
614 // Tracks are fitted only within
615 // the range defined by the user.
616 Int_t nModules = 0;
98f0f87b 617 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
25be1e5c 618 nModules += AliGeomManager::LayerSize(iLayer);
cc345ce3 619 TArrayI volIds(nModules);
620
621 Int_t modnum = 0;
98f0f87b 622 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
25be1e5c 623 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
624 UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
cc345ce3 625 volIds.AddAt(volId,modnum);
626 modnum++;
627 }
628 }
629
98937d93 630 while (iterations > 0) {
cc345ce3 631 AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
632 iterations--;
98937d93 633 }
634}
635
636//______________________________________________________________________________
25be1e5c 637void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
638 AliGeomManager::ELayerID layerRangeMin,
639 AliGeomManager::ELayerID layerRangeMax,
98937d93 640 Int_t iterations)
641{
642 // Align detector volumes within
643 // a given layer.
644 // Tracks are fitted only within
645 // the range defined by the user.
25be1e5c 646 Int_t nModules = AliGeomManager::LayerSize(layer);
cc345ce3 647 TArrayI volIds(nModules);
648 for (Int_t iModule = 0; iModule < nModules; iModule++) {
25be1e5c 649 UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
cc345ce3 650 volIds.AddAt(volId,iModule);
651 }
652
98937d93 653 while (iterations > 0) {
cc345ce3 654 AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
655 iterations--;
656 }
657}
658
659//______________________________________________________________________________
660void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
661 Int_t iterations)
662{
663 // Align single detector volume to
664 // another volume.
665 // Tracks are fitted only within
666 // the second volume.
667 TArrayI volIds(1);
668 volIds.AddAt(volId,0);
669 TArrayI volIdsFit(1);
670 volIdsFit.AddAt(volIdFit,0);
671
672 while (iterations > 0) {
673 AlignVolumes(&volIds,&volIdsFit);
98937d93 674 iterations--;
675 }
676}
677
678//______________________________________________________________________________
cc345ce3 679void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
25be1e5c 680 AliGeomManager::ELayerID layerRangeMin,
681 AliGeomManager::ELayerID layerRangeMax,
46ae650f 682 Int_t iterations)
98937d93 683{
cc345ce3 684 // Align a set of detector volumes.
98937d93 685 // Tracks are fitted only within
46ae650f 686 // the range defined by the user
687 // (by layerRangeMin and layerRangeMax)
cc345ce3 688 // or within the set of volidsfit
46ae650f 689 // Repeat the procedure 'iterations' times
98937d93 690
cc345ce3 691 Int_t nVolIds = volids->GetSize();
692 if (nVolIds == 0) {
693 AliError("Volume IDs array is empty!");
694 return;
695 }
98937d93 696
cc345ce3 697 // Load only the tracks with at least one
698 // space point in the set of volume (volids)
98937d93 699 BuildIndex();
700 AliTrackPointArray **points;
46ae650f 701 // Start the iterations
702 while (iterations > 0) {
cc345ce3 703 Int_t nArrays = LoadPoints(volids, points);
46ae650f 704 if (nArrays == 0) return;
705
706 AliTrackResiduals *minimizer = CreateMinimizer();
707 minimizer->SetNTracks(nArrays);
cc345ce3 708 minimizer->InitAlignObj();
46ae650f 709 AliTrackFitter *fitter = CreateFitter();
710 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
33bebad6 711 if (!points[iArray]) continue;
46ae650f 712 fitter->SetTrackPointArray(points[iArray], kFALSE);
33bebad6 713 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
46ae650f 714 AliTrackPointArray *pVolId,*pTrack;
715 fitter->GetTrackResiduals(pVolId,pTrack);
716 minimizer->AddTrackPointArrays(pVolId,pTrack);
717 }
718 minimizer->Minimize();
cc345ce3 719
720 // Update the alignment object(s)
91264393 721 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
cc345ce3 722 UShort_t volid = (*volids)[iVolId];
723 Int_t iModule;
25be1e5c 724 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
725 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 726 *alignObj *= *minimizer->GetAlignObj();
727 alignObj->Print("");
728 }
98937d93 729
46ae650f 730 UnloadPoints(nArrays, points);
731
732 iterations--;
98937d93 733 }
98937d93 734}
735
736//______________________________________________________________________________
cc345ce3 737Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
98937d93 738{
739 // Load track point arrays with at least
cc345ce3 740 // one space point in a given set of detector
741 // volumes (array volids).
98937d93 742 // Use the already created tree index for
743 // fast access.
98937d93 744
46ae650f 745 if (!fPointsTree) {
cc345ce3 746 AliError("Tree with the space point arrays not initialized!");
98937d93 747 points = 0;
748 return 0;
749 }
750
cc345ce3 751 Int_t nVolIds = volids->GetSize();
752 if (nVolIds == 0) {
753 AliError("Volume IDs array is empty!");
754 points = 0;
755 return 0;
756 }
757
758 Int_t nArrays = 0;
759 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
760 UShort_t volid = (*volids)[iVolId];
761 Int_t iModule;
25be1e5c 762 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
cc345ce3 763
764 // In case of empty index
25be1e5c 765 if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
cc345ce3 766 AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
767 continue;
768 }
25be1e5c 769 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 770 }
46ae650f 771
46ae650f 772 if (nArrays == 0) {
cc345ce3 773 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
98937d93 774 points = 0;
775 return 0;
776 }
777
98937d93 778 AliTrackPointArray* array = 0;
779 fPointsTree->SetBranchAddress("SP", &array);
780
cc345ce3 781 // Allocate the pointer to the space-point arrays
98937d93 782 points = new AliTrackPointArray*[nArrays];
5caa5567 783 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
cc345ce3 784
785 // Init the array used to flag already loaded tree entries
7fa594d2 786 Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
cc345ce3 787 for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
788 indexUsed[i] = kFALSE;
789
790 // Start the loop over the volume ids
791 Int_t iArray = 0;
792 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
793 UShort_t volid = (*volids)[iVolId];
794 Int_t iModule;
25be1e5c 795 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
cc345ce3 796
25be1e5c 797 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
798 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 799 AliTrackPoint p;
800
801 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
802
803 // Get tree entry
804 Int_t entry = (*index)[iArrayId];
805 if (indexUsed[entry] == kTRUE) continue;
806 fPointsTree->GetEvent(entry);
807 if (!array) {
808 AliWarning("Wrong space point array index!");
809 continue;
46ae650f 810 }
cc345ce3 811 indexUsed[entry] = kTRUE;
812
813 // Get the space-point array
814 Int_t nPoints = array->GetNPoints();
815 points[iArray] = new AliTrackPointArray(nPoints);
816 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
817 array->GetPoint(p,iPoint);
818 Int_t modnum;
25be1e5c 819 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
33bebad6 820 // check if the layer id is valid
25be1e5c 821 if ((layer < AliGeomManager::kFirstLayer) ||
822 (layer >= AliGeomManager::kLastLayer)) {
6baf24ce 823 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
25be1e5c 824 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
6baf24ce 825 continue;
826 }
25be1e5c 827 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
6baf24ce 828 (modnum < 0)) {
829 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
25be1e5c 830 layer,modnum,AliGeomManager::LayerSize(layer)));
6baf24ce 831 continue;
832 }
cc345ce3 833
834 // Misalignment is introduced here
835 // Switch it off in case of real
836 // alignment job!
837 if (fMisalignObjs) {
25be1e5c 838 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
cc345ce3 839 if (misalignObj)
840 misalignObj->Transform(p);
841 }
842 // End of misalignment
46ae650f 843
25be1e5c 844 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
cc345ce3 845 alignObj->Transform(p);
846 points[iArray]->AddPoint(iPoint,&p);
847 }
848 iArray++;
98937d93 849 }
850 }
851
cc345ce3 852
853 delete [] indexUsed;
854
98937d93 855 return nArrays;
856}
857
858//______________________________________________________________________________
859void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
860{
861 // Unload track point arrays for a given
862 // detector volume
863 for (Int_t iArray = 0; iArray < n; iArray++)
864 delete points[iArray];
865 delete [] points;
866}
867
868//______________________________________________________________________________
869AliTrackFitter *AliAlignmentTracks::CreateFitter()
870{
871 // Check if the user has already supplied
872 // a track fitter object.
873 // If not, create a default one.
874 if (!fTrackFitter)
875 fTrackFitter = new AliTrackFitterRieman;
876
877 return fTrackFitter;
878}
879
880//______________________________________________________________________________
881AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
882{
883 // Check if the user has already supplied
884 // a track residuals minimizer object.
885 // If not, create a default one.
886 if (!fMinimizer)
887 fMinimizer = new AliTrackResidualsChi2;
888
889 return fMinimizer;
890}
46ae650f 891
892//______________________________________________________________________________
893Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
894{
895 // The method reads from a file a set of AliAlignObj which are
896 // then used to apply misalignments directly on the track
897 // space-points. The method is supposed to be used only for
898 // fast development and debugging of the alignment algorithms.
899 // Be careful not to use it in the case of 'real' alignment
900 // scenario since it will bias the results.
901
902 // Initialize the misalignment objects array
25be1e5c 903 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
46ae650f 904 fMisalignObjs = new AliAlignObj**[nLayers];
25be1e5c 905 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
906 fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
907 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
46ae650f 908 fMisalignObjs[iLayer][iModule] = 0x0;
909 }
910
911 // Open the misliagnment file and load the array with
912 // misalignment objects
913 TFile* inFile = TFile::Open(misalignObjFileName,"READ");
914 if (!inFile || !inFile->IsOpen()) {
915 AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
916 return kFALSE;
917 }
918
919 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
920 if (!array) {
921 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
922 inFile->Close();
923 return kFALSE;
924 }
925 inFile->Close();
926
927 // Store the misalignment objects for further usage
928 Int_t nObjs = array->GetEntriesFast();
25be1e5c 929 AliGeomManager::ELayerID layerId; // volume layer
46ae650f 930 Int_t modId; // volume ID inside the layer
931 for(Int_t i=0; i<nObjs; i++)
932 {
933 AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
934 alObj->GetVolUID(layerId,modId);
25be1e5c 935 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
46ae650f 936 }
937 return kTRUE;
938}