New plots in QA plotting macro
[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"
4fc70c8b 36#include "AliESDfriend.h"
98937d93 37#include "AliLog.h"
98937d93 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
465//______________________________________________________________________________
98937d93 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);
6621b37d 775 if (nArrays == 0) {
776 UnloadPoints(pointsdim, points);
777 return kFALSE;
778 }
46ae650f 779
780 AliTrackResiduals *minimizer = CreateMinimizer();
781 minimizer->SetNTracks(nArrays);
cc345ce3 782 minimizer->InitAlignObj();
46ae650f 783 AliTrackFitter *fitter = CreateFitter();
784 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
33bebad6 785 if (!points[iArray]) continue;
46ae650f 786 fitter->SetTrackPointArray(points[iArray], kFALSE);
33bebad6 787 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
46ae650f 788 AliTrackPointArray *pVolId,*pTrack;
789 fitter->GetTrackResiduals(pVolId,pTrack);
790 minimizer->AddTrackPointArrays(pVolId,pTrack);
791 }
5ffa0b18 792 if (!(result = minimizer->Minimize())) {
793 UnloadPoints(pointsdim, points);
794 break;
795 }
cc345ce3 796
797 // Update the alignment object(s)
91264393 798 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
cc345ce3 799 UShort_t volid = (*volids)[iVolId];
800 Int_t iModule;
25be1e5c 801 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
802 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 803 *alignObj *= *minimizer->GetAlignObj();
c2c53eb8 804 if(iterations==1)alignObj->Print("");
cc345ce3 805 }
98937d93 806
a5b84e94 807 UnloadPoints(pointsdim, points);
46ae650f 808
809 iterations--;
98937d93 810 }
6fdedff6 811 return result;
98937d93 812}
813
814//______________________________________________________________________________
a5b84e94 815Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points,Int_t &pointsdim)
98937d93 816{
817 // Load track point arrays with at least
cc345ce3 818 // one space point in a given set of detector
819 // volumes (array volids).
98937d93 820 // Use the already created tree index for
821 // fast access.
98937d93 822
46ae650f 823 if (!fPointsTree) {
cc345ce3 824 AliError("Tree with the space point arrays not initialized!");
98937d93 825 points = 0;
826 return 0;
827 }
828
cc345ce3 829 Int_t nVolIds = volids->GetSize();
830 if (nVolIds == 0) {
831 AliError("Volume IDs array is empty!");
832 points = 0;
833 return 0;
834 }
835
836 Int_t nArrays = 0;
837 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
838 UShort_t volid = (*volids)[iVolId];
839 Int_t iModule;
25be1e5c 840 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
cc345ce3 841
842 // In case of empty index
25be1e5c 843 if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
cc345ce3 844 AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
845 continue;
846 }
25be1e5c 847 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 848 }
46ae650f 849
46ae650f 850 if (nArrays == 0) {
cc345ce3 851 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
6b7c7eba 852 points = 0x0;
98937d93 853 return 0;
854 }
855
98937d93 856 AliTrackPointArray* array = 0;
857 fPointsTree->SetBranchAddress("SP", &array);
858
cc345ce3 859 // Allocate the pointer to the space-point arrays
a5b84e94 860 pointsdim=nArrays;
98937d93 861 points = new AliTrackPointArray*[nArrays];
5caa5567 862 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
cc345ce3 863
864 // Init the array used to flag already loaded tree entries
7fa594d2 865 Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
cc345ce3 866 for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
867 indexUsed[i] = kFALSE;
868
869 // Start the loop over the volume ids
870 Int_t iArray = 0;
871 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
872 UShort_t volid = (*volids)[iVolId];
873 Int_t iModule;
25be1e5c 874 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
cc345ce3 875
25be1e5c 876 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
877 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
cc345ce3 878 AliTrackPoint p;
879
880 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
881
882 // Get tree entry
883 Int_t entry = (*index)[iArrayId];
8069080f 884 if (indexUsed[entry] == kTRUE) {
885 nArrays--;
886 continue;
887 }
cc345ce3 888 fPointsTree->GetEvent(entry);
889 if (!array) {
890 AliWarning("Wrong space point array index!");
891 continue;
46ae650f 892 }
cc345ce3 893 indexUsed[entry] = kTRUE;
894
895 // Get the space-point array
896 Int_t nPoints = array->GetNPoints();
897 points[iArray] = new AliTrackPointArray(nPoints);
898 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
899 array->GetPoint(p,iPoint);
900 Int_t modnum;
25be1e5c 901 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
33bebad6 902 // check if the layer id is valid
25be1e5c 903 if ((layer < AliGeomManager::kFirstLayer) ||
904 (layer >= AliGeomManager::kLastLayer)) {
6baf24ce 905 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
25be1e5c 906 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
6baf24ce 907 continue;
908 }
25be1e5c 909 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
6baf24ce 910 (modnum < 0)) {
911 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
25be1e5c 912 layer,modnum,AliGeomManager::LayerSize(layer)));
6baf24ce 913 continue;
914 }
cc345ce3 915
916 // Misalignment is introduced here
917 // Switch it off in case of real
918 // alignment job!
919 if (fMisalignObjs) {
25be1e5c 920 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
cc345ce3 921 if (misalignObj)
922 misalignObj->Transform(p);
923 }
924 // End of misalignment
46ae650f 925
6b7c7eba 926
25be1e5c 927 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
6b7c7eba 928 UShort_t volp=p.GetVolumeID();
929 Bool_t found=kFALSE;
930 if(fCovIsUsed){
931 for (Int_t iVol = 0; iVol < nVolIds; iVol++) {
932 UShort_t vol = (*volids)[iVol];
933 if(volp==vol){
934 alignObj->Transform(p,kFALSE);
935 found=kTRUE;
936 break;
937 }
938 }
939 }
940 if(!found)alignObj->Transform(p,fCovIsUsed);
cc345ce3 941 points[iArray]->AddPoint(iPoint,&p);
942 }
943 iArray++;
98937d93 944 }
945 }
946
cc345ce3 947
948 delete [] indexUsed;
949
98937d93 950 return nArrays;
951}
952
953//______________________________________________________________________________
954void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
955{
956 // Unload track point arrays for a given
957 // detector volume
958 for (Int_t iArray = 0; iArray < n; iArray++)
959 delete points[iArray];
960 delete [] points;
961}
962
963//______________________________________________________________________________
964AliTrackFitter *AliAlignmentTracks::CreateFitter()
965{
966 // Check if the user has already supplied
967 // a track fitter object.
968 // If not, create a default one.
969 if (!fTrackFitter)
970 fTrackFitter = new AliTrackFitterRieman;
971
972 return fTrackFitter;
973}
974
975//______________________________________________________________________________
976AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
977{
978 // Check if the user has already supplied
979 // a track residuals minimizer object.
980 // If not, create a default one.
981 if (!fMinimizer)
982 fMinimizer = new AliTrackResidualsChi2;
983
984 return fMinimizer;
985}
46ae650f 986
987//______________________________________________________________________________
988Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
989{
990 // The method reads from a file a set of AliAlignObj which are
991 // then used to apply misalignments directly on the track
992 // space-points. The method is supposed to be used only for
993 // fast development and debugging of the alignment algorithms.
994 // Be careful not to use it in the case of 'real' alignment
995 // scenario since it will bias the results.
996
997 // Initialize the misalignment objects array
25be1e5c 998 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
46ae650f 999 fMisalignObjs = new AliAlignObj**[nLayers];
25be1e5c 1000 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
1001 fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
1002 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
46ae650f 1003 fMisalignObjs[iLayer][iModule] = 0x0;
1004 }
1005
1006 // Open the misliagnment file and load the array with
1007 // misalignment objects
1008 TFile* inFile = TFile::Open(misalignObjFileName,"READ");
1009 if (!inFile || !inFile->IsOpen()) {
1010 AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
1011 return kFALSE;
1012 }
1013
1014 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
1015 if (!array) {
1016 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
1017 inFile->Close();
1018 return kFALSE;
1019 }
1020 inFile->Close();
1021
1022 // Store the misalignment objects for further usage
1023 Int_t nObjs = array->GetEntriesFast();
25be1e5c 1024 AliGeomManager::ELayerID layerId; // volume layer
46ae650f 1025 Int_t modId; // volume ID inside the layer
1026 for(Int_t i=0; i<nObjs; i++)
1027 {
1028 AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
1029 alObj->GetVolUID(layerId,modId);
6b7c7eba 1030 if(layerId<AliGeomManager::kFirstLayer) {
1031 AliWarning(Form("Alignment object is ignored: %s",alObj->GetSymName()));
1032 continue;
1033 }
25be1e5c 1034 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
46ae650f 1035 }
1036 return kTRUE;
1037}
d64e3399 1038
1039
1040//________________________________________________
1041void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
1042
1043 Int_t last=0;
1044 TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
1045 TClonesArray &alo=*clonesArray;
1046 for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
1047
1048 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
1049
1050 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
1051 new(alo[last])AliAlignObjParams(*alignObj);
1052 last++;
1053 }
1054 }
1055 TFile *file=new TFile(outfilename.Data(),"RECREATE");
1056 file->cd();
1057
1058 alo.Write("ITSAlignObjs",TObject::kSingleKey);
1059 file->Close();
1060
1061
1062 delete clonesArray;
1063 return;
1064}