]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - STEER/AliAlignmentTracks.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 steering class ro run QA stand alone
[u/mrichter/AliRoot.git] / STEER / AliAlignmentTracks.cxx
... / ...
CommitLineData
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>
26#include <TFile.h>
27#include <TVector3.h>
28
29#include "AliAlignmentTracks.h"
30#include "AliTrackPointArray.h"
31#include "AliAlignObjParams.h"
32#include "AliTrackFitterRieman.h"
33#include "AliTrackResidualsChi2.h"
34#include "AliESDEvent.h"
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),
48 fAlignObjs(0),
49 fMisalignObjs(0),
50 fTrackFitter(0),
51 fMinimizer(0),
52 fDoUpdate(kTRUE)
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),
68 fAlignObjs(0),
69 fMisalignObjs(0),
70 fTrackFitter(0),
71 fMinimizer(0),
72 fDoUpdate(kTRUE)
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):
84 fESDChain(new TChain(esdtreename)),
85 fPointsFilename("AliTrackPoints.root"),
86 fPointsFile(0),
87 fPointsTree(0),
88 fLastIndex(0),
89 fArrayIndex(0),
90 fIsIndexBuilt(kFALSE),
91 fAlignObjs(0),
92 fMisalignObjs(0),
93 fTrackFitter(0),
94 fMinimizer(0),
95 fDoUpdate(kTRUE)
96{
97 // Constructor in the case
98 // the user provides a single ESD file
99 // or a directory containing ESD files
100 fESDChain->Add(esdfilename);
101
102 InitIndex();
103 InitAlignObjs();
104}
105
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
146//________________________________________________________________________
147void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
148 Int_t minITSpts,
149 Bool_t cuts,
150 Float_t minAngleWrtITSModulePlanes,
151 Float_t minMom,Float_t maxMom,
152 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
153 Float_t minSinTheta,Float_t maxSinTheta)
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
161 AliESDEvent *esd = new AliESDEvent();
162 esd->ReadFromTree(fESDChain);
163 AliESDfriend *esdf = 0;
164 fESDChain->SetBranchStatus("ESDfriend*",1);
165 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
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");
180 const AliTrackPointArray *array = 0;
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
188 Int_t ievent = 0;
189 while (fESDChain->GetEntry(ievent++)) {
190 if (!esd) break;
191
192 esd->SetESDfriend(esdf); //Attach the friend to the ESD
193
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;
198
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 }
207
208 AliTrackPoint point;
209 array = track->GetTrackPointArray();
210
211 if(onlyITS) {
212 array2 = new AliTrackPointArray(track->GetNcls(0));
213 Bool_t smallAngleWrtITSModulePlanes=kFALSE;
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 }
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 }
235 }
236 if(smallAngleWrtITSModulePlanes) continue;
237 } // end if(onlyITS)
238
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();
249
250 return;
251}
252
253//_____________________________________________________________________________
254void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
255 Int_t minITSpts,Float_t maxMatchingAngle,
256 Bool_t cuts,
257 Float_t minAngleWrtITSModulePlanes,
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
269 AliESDEvent *esd = new AliESDEvent();
270 esd->ReadFromTree(fESDChain);
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();
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];
303 Int_t ngt=0;
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;
309 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
310 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
311 if(cuts) {
312 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
313 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
314 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
315 Float_t sintheta = TMath::Sin(theta);
316 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
317 }
318 goodtracksArray[ngt]=itrack;
319 phiArray[ngt]=phi;
320 thetaArray[ngt]=theta;
321 ngt++;
322 }
323
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;
340 //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
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;
354
355 AliESDtrack * track1 = esd->GetTrack(good1);
356 AliESDtrack * track2 = esd->GetTrack(good2);
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;
369 Bool_t smallAngleWrtITSModulePlanes=kFALSE;
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;
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 }
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;
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 }
410 }
411 array2->AddPoint(jpt,&point);
412 jpt++;
413 }
414
415 if(smallAngleWrtITSModulePlanes) continue;
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();
425 return;
426}
427
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
439 if (fIsIndexBuilt) return;
440
441 fIsIndexBuilt = kTRUE;
442
443 // Dummy object is created in order
444 // to initialize the volume paths
445 AliAlignObjParams alobj;
446
447 fPointsFile = TFile::Open(fPointsFilename);
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;
455 fPointsTree = (TTree*) fPointsFile->Get("spTree");
456 if (!fPointsTree) {
457 AliWarning("No pointsTree found!");
458 return;
459 }
460 fPointsTree->SetBranchAddress("SP", &array);
461
462 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
463 for (Int_t iArray = 0; iArray < nArrays; iArray++)
464 {
465 fPointsTree->GetEvent(iArray);
466 if (!array) continue;
467 for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
468 UShort_t volId = array->GetVolumeID()[ipoint];
469 // check if the volId is valid
470 if (!AliGeomManager::SymName(volId)) {
471 AliError(Form("The volume id %d has no default volume name !",
472 volId));
473 continue;
474 }
475 Int_t modId;
476 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
477 - AliGeomManager::kFirstLayer;
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
505//______________________________________________________________________________
506void AliAlignmentTracks::InitIndex()
507{
508 // Initialize the index arrays
509 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
510 fLastIndex = new Int_t*[nLayers];
511 fArrayIndex = new TArrayI**[nLayers];
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++) {
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
527
528 fIsIndexBuilt = kFALSE;
529
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++) {
532 fLastIndex[iLayer][iModule] = 0;
533 }
534 }
535}
536
537//______________________________________________________________________________
538void AliAlignmentTracks::DeleteIndex()
539{
540 // Delete the index arrays
541 // Called by the destructor
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++) {
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//______________________________________________________________________________
557Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
558{
559 // Read alignment object from a file
560 // To be replaced by a call to CDB
561 AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
562
563 return kFALSE;
564}
565
566//______________________________________________________________________________
567void AliAlignmentTracks::InitAlignObjs()
568{
569 // Initialize the alignment objects array
570 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
571 fAlignObjs = new AliAlignObj**[nLayers];
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);
576 fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
577 }
578 }
579}
580
581//______________________________________________________________________________
582void AliAlignmentTracks::ResetAlignObjs()
583{
584 // Reset the alignment objects array
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++)
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
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++)
597 if (fAlignObjs[iLayer][iModule])
598 delete fAlignObjs[iLayer][iModule];
599 delete [] fAlignObjs[iLayer];
600 }
601 delete [] fAlignObjs;
602 fAlignObjs = 0;
603}
604
605void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
606 AliGeomManager::ELayerID lastLayer,
607 AliGeomManager::ELayerID layerRangeMin,
608 AliGeomManager::ELayerID layerRangeMax,
609 Int_t iterations)
610{
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;
617 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
618 nModules += AliGeomManager::LayerSize(iLayer);
619 TArrayI volIds(nModules);
620
621 Int_t modnum = 0;
622 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
623 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
624 UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
625 volIds.AddAt(volId,modnum);
626 modnum++;
627 }
628 }
629
630 while (iterations > 0) {
631 AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
632 iterations--;
633 }
634}
635
636//______________________________________________________________________________
637void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
638 AliGeomManager::ELayerID layerRangeMin,
639 AliGeomManager::ELayerID layerRangeMax,
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.
646 Int_t nModules = AliGeomManager::LayerSize(layer);
647 TArrayI volIds(nModules);
648 for (Int_t iModule = 0; iModule < nModules; iModule++) {
649 UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
650 volIds.AddAt(volId,iModule);
651 }
652
653 while (iterations > 0) {
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);
674 iterations--;
675 }
676}
677
678//______________________________________________________________________________
679void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
680 AliGeomManager::ELayerID layerRangeMin,
681 AliGeomManager::ELayerID layerRangeMax,
682 Int_t iterations)
683{
684 // Align a set of detector volumes.
685 // Tracks are fitted only within
686 // the range defined by the user
687 // (by layerRangeMin and layerRangeMax)
688 // or within the set of volidsfit
689 // Repeat the procedure 'iterations' times
690
691 Int_t nVolIds = volids->GetSize();
692 if (nVolIds == 0) {
693 AliError("Volume IDs array is empty!");
694 return;
695 }
696
697 // Load only the tracks with at least one
698 // space point in the set of volume (volids)
699 BuildIndex();
700 AliTrackPointArray **points;
701 // Start the iterations
702 while (iterations > 0) {
703 Int_t nArrays = LoadPoints(volids, points);
704 if (nArrays == 0) return;
705
706 AliTrackResiduals *minimizer = CreateMinimizer();
707 minimizer->SetNTracks(nArrays);
708 minimizer->InitAlignObj();
709 AliTrackFitter *fitter = CreateFitter();
710 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
711 if (!points[iArray]) continue;
712 fitter->SetTrackPointArray(points[iArray], kFALSE);
713 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
714 AliTrackPointArray *pVolId,*pTrack;
715 fitter->GetTrackResiduals(pVolId,pTrack);
716 minimizer->AddTrackPointArrays(pVolId,pTrack);
717 }
718 minimizer->Minimize();
719
720 // Update the alignment object(s)
721 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
722 UShort_t volid = (*volids)[iVolId];
723 Int_t iModule;
724 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
725 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
726 *alignObj *= *minimizer->GetAlignObj();
727 alignObj->Print("");
728 }
729
730 UnloadPoints(nArrays, points);
731
732 iterations--;
733 }
734}
735
736//______________________________________________________________________________
737Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
738{
739 // Load track point arrays with at least
740 // one space point in a given set of detector
741 // volumes (array volids).
742 // Use the already created tree index for
743 // fast access.
744
745 if (!fPointsTree) {
746 AliError("Tree with the space point arrays not initialized!");
747 points = 0;
748 return 0;
749 }
750
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;
762 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
763
764 // In case of empty index
765 if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
766 AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
767 continue;
768 }
769 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
770 }
771
772 if (nArrays == 0) {
773 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
774 points = 0;
775 return 0;
776 }
777
778 AliTrackPointArray* array = 0;
779 fPointsTree->SetBranchAddress("SP", &array);
780
781 // Allocate the pointer to the space-point arrays
782 points = new AliTrackPointArray*[nArrays];
783 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
784
785 // Init the array used to flag already loaded tree entries
786 Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
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;
795 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
796
797 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
798 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
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;
810 }
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;
819 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
820 // check if the layer id is valid
821 if ((layer < AliGeomManager::kFirstLayer) ||
822 (layer >= AliGeomManager::kLastLayer)) {
823 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
824 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
825 continue;
826 }
827 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
828 (modnum < 0)) {
829 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
830 layer,modnum,AliGeomManager::LayerSize(layer)));
831 continue;
832 }
833
834 // Misalignment is introduced here
835 // Switch it off in case of real
836 // alignment job!
837 if (fMisalignObjs) {
838 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
839 if (misalignObj)
840 misalignObj->Transform(p);
841 }
842 // End of misalignment
843
844 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
845 alignObj->Transform(p);
846 points[iArray]->AddPoint(iPoint,&p);
847 }
848 iArray++;
849 }
850 }
851
852
853 delete [] indexUsed;
854
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}
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
903 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
904 fMisalignObjs = new AliAlignObj**[nLayers];
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++)
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();
929 AliGeomManager::ELayerID layerId; // volume layer
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);
935 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
936 }
937 return kTRUE;
938}