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