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