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