Alignment framework (C.Cheshkov). More information is available in http://agenda...
[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>
26
27#include "AliAlignmentTracks.h"
28#include "AliTrackPointArray.h"
29#include "AliAlignObjAngles.h"
30#include "AliTrackFitterRieman.h"
31#include "AliTrackResidualsChi2.h"
32#include "AliESD.h"
33#include "AliLog.h"
34
35ClassImp(AliAlignmentTracks)
36
37//______________________________________________________________________________
38AliAlignmentTracks::AliAlignmentTracks():
39 fESDChain(0),
40 fPointsFilename("AliTrackPoints.root"),
41 fPointsFile(0),
42 fPointsTree(0),
43 fLastIndex(0),
44 fArrayIndex(0),
45 fIsIndexBuilt(kFALSE),
46 fTrackFitter(0),
47 fMinimizer(0)
48{
49 // Default constructor
50 InitIndex();
51 InitAlignObjs();
52}
53
54//______________________________________________________________________________
55AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
56 fESDChain(esdchain),
57 fPointsFilename("AliTrackPoints.root"),
58 fPointsFile(0),
59 fPointsTree(0),
60 fLastIndex(0),
61 fArrayIndex(0),
62 fIsIndexBuilt(kFALSE),
63 fTrackFitter(0),
64 fMinimizer(0)
65{
66 // Constructor in the case
67 // the user provides an already
68 // built TChain with ESD trees
69 InitIndex();
70 InitAlignObjs();
71}
72
73
74//______________________________________________________________________________
75AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
76 fPointsFilename("AliTrackPoints.root"),
77 fPointsFile(0),
78 fPointsTree(0),
79 fLastIndex(0),
80 fArrayIndex(0),
81 fIsIndexBuilt(kFALSE),
82 fTrackFitter(0),
83 fMinimizer(0)
84{
85 // Constructor in the case
86 // the user provides a single ESD file
87 // or a directory containing ESD files
88 fESDChain = new TChain(esdtreename);
89 fESDChain->Add(esdfilename);
90
91 InitIndex();
92 InitAlignObjs();
93}
94
95//______________________________________________________________________________
96AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
97 TObject(alignment)
98{
99 // Copy constructor
100 // not implemented
101 AliWarning("Copy constructor not implemented!");
102}
103
104//______________________________________________________________________________
105AliAlignmentTracks& AliAlignmentTracks::operator= (const AliAlignmentTracks& alignment)
106{
107 // Asignment operator
108 // not implemented
109 if(this==&alignment) return *this;
110
111 AliWarning("Asignment operator not implemented!");
112
113 ((TObject *)this)->operator=(alignment);
114
115 return *this;
116}
117
118//______________________________________________________________________________
119AliAlignmentTracks::~AliAlignmentTracks()
120{
121 // Destructor
122 if (fESDChain) delete fESDChain;
123
124 DeleteIndex();
125 DeleteAlignObjs();
126
127 delete fTrackFitter;
128 delete fMinimizer;
129
130 if (fPointsFile) fPointsFile->Close();
131}
132
133//______________________________________________________________________________
134void AliAlignmentTracks::AddESD(TChain *esdchain)
135{
136 // Add a chain with ESD files
137 if (fESDChain)
138 fESDChain->Add(esdchain);
139 else
140 fESDChain = esdchain;
141}
142
143//______________________________________________________________________________
144void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
145{
146 // Add a single file or
147 // a directory to the chain
148 // with the ESD files
149 if (fESDChain)
150 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
151 else {
152 fESDChain = new TChain(esdtreename);
153 fESDChain->Add(esdfilename);
154 }
155}
156
157//______________________________________________________________________________
158void AliAlignmentTracks::ProcessESD()
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
166 AliESD *esd = 0;
167 fESDChain->SetBranchAddress("ESD",&esd);
168
169 // Open the output file
170 if (fPointsFilename.Data() == "") {
171 AliWarning("Incorrect output filename!");
172 return;
173 }
174
175 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
176 if (!pointsFile || !pointsFile->IsOpen()) {
177 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
178 return;
179 }
180
181 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
182 AliTrackPointArray *array = 0;
183 pointsTree->Branch("SP","AliTrackPointArray", &array);
184
185 Int_t ievent = 0;
186 while (fESDChain->GetEntry(ievent++)) {
187 if (!esd) break;
188 Int_t ntracks = esd->GetNumberOfTracks();
189 for (Int_t itrack=0; itrack < ntracks; itrack++) {
190 AliESDtrack * track = esd->GetTrack(itrack);
191 if (!track) continue;
192
193 UInt_t status = AliESDtrack::kITSpid;
194 status|=AliESDtrack::kTPCpid;
195 status|=AliESDtrack::kTRDpid;
196 if ((track->GetStatus() & status) != status) continue;
197
198 if (track->GetP() < 0.5) continue;
199
200 array = track->GetTrackPointArray();
201 pointsTree->Fill();
202 }
203 }
204
205 if (!pointsTree->Write()) {
206 AliWarning("Can't write the tree with track point arrays!");
207 return;
208 }
209
210 pointsFile->Close();
211}
212
213//______________________________________________________________________________
214void AliAlignmentTracks::ProcessESD(TSelector *selector)
215{
216 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
217}
218
219//______________________________________________________________________________
220void AliAlignmentTracks::BuildIndex()
221{
222 // Build index of points tree entries
223 // Used for access based on the volume IDs
224 if (fIsIndexBuilt)
225 ResetIndex();
226 else
227 fIsIndexBuilt = kTRUE;
228
229 TFile *fPointsFile = TFile::Open(fPointsFilename);
230 if (!fPointsFile || !fPointsFile->IsOpen()) {
231 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
232 return;
233 }
234
235 // AliTrackPointArray* array = new AliTrackPointArray;
236 AliTrackPointArray* array = 0;
237 TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
238 if (!pointsTree) {
239 AliWarning("No pointsTree found!");
240 return;
241 }
242 pointsTree->SetBranchAddress("SP", &array);
243
244 Int_t nArrays = pointsTree->GetEntries();
245 for (Int_t iArray = 0; iArray < nArrays; iArray++)
246 {
247 pointsTree->GetEvent(iArray);
248 if (!array) continue;
249 for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
250 UShort_t volId = array->GetVolumeID()[ipoint];
251 Int_t modId;
252 Int_t layerId = AliAlignObj::VolUIDToLayer(volId,modId)
253 - AliAlignObj::kFirstLayer;
254 if (!fArrayIndex[layerId][modId]) {
255 //first entry for this volume
256 fArrayIndex[layerId][modId] = new TArrayI(1000);
257 }
258 else {
259 Int_t size = fArrayIndex[layerId][modId]->GetSize();
260 // If needed allocate new size
261 if (fLastIndex[layerId][modId] >= size)
262 fArrayIndex[layerId][modId]->Set(size + 1000);
263 }
264
265 // Check if the index is already filled
266 Bool_t fillIndex = kTRUE;
267 if (fLastIndex[layerId][modId] != 0) {
268 if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
269 fillIndex = kFALSE;
270 }
271 // Fill the index array and store last filled index
272 if (fillIndex) {
273 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
274 fLastIndex[layerId][modId]++;
275 }
276 }
277 }
278
279}
280
281// //______________________________________________________________________________
282// void AliAlignmentTracks::BuildIndexLayer(AliAlignObj::ELayerID layer)
283// {
284// }
285
286// //______________________________________________________________________________
287// void AliAlignmentTracks::BuildIndexVolume(UShort_t volid)
288// {
289// }
290
291//______________________________________________________________________________
292void AliAlignmentTracks::InitIndex()
293{
294 // Initialize the index arrays
295 Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
296 fLastIndex = new Int_t*[nLayers];
297 fArrayIndex = new TArrayI**[nLayers];
298 for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
299 fLastIndex[iLayer] = new Int_t[AliAlignObj::LayerSize(iLayer)];
300 fArrayIndex[iLayer] = new TArrayI*[AliAlignObj::LayerSize(iLayer)];
301 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
302 fLastIndex[iLayer][iModule] = 0;
303 fArrayIndex[iLayer][iModule] = 0;
304 }
305 }
306}
307
308//______________________________________________________________________________
309void AliAlignmentTracks::ResetIndex()
310{
311 // Reset the value of the last filled index
312 // Do not realocate memory
313 for (Int_t iLayer = AliAlignObj::kFirstLayer; iLayer < AliAlignObj::kLastLayer; iLayer++) {
314 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
315 fLastIndex[iLayer][iModule] = 0;
316 }
317 }
318}
319
320//______________________________________________________________________________
321void AliAlignmentTracks::DeleteIndex()
322{
323 // Delete the index arrays
324 // Called by the destructor
325 for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
326 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
327 if (fArrayIndex[iLayer][iModule]) {
328 delete fArrayIndex[iLayer][iModule];
329 fArrayIndex[iLayer][iModule] = 0;
330 }
331 }
332 delete [] fLastIndex[iLayer];
333 delete [] fArrayIndex[iLayer];
334 }
335 delete [] fLastIndex;
336 delete [] fArrayIndex;
337}
338
339//______________________________________________________________________________
340Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
341{
342 // Read alignment object from a file
343 // To be replaced by a call to CDB
344 AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));
345
346 return kFALSE;
347}
348
349//______________________________________________________________________________
350void AliAlignmentTracks::InitAlignObjs()
351{
352 // Initialize the alignment objects array
353 Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
354 fAlignObjs = new AliAlignObj**[nLayers];
355 for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
356 fAlignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
357 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
358 fAlignObjs[iLayer][iModule] = new AliAlignObjAngles;
359 }
360}
361
362//______________________________________________________________________________
363void AliAlignmentTracks::ResetAlignObjs()
364{
365 // Reset the alignment objects array
366 for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
367 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
368 fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
369 }
370}
371
372//______________________________________________________________________________
373void AliAlignmentTracks::DeleteAlignObjs()
374{
375 // Delete the alignment objects array
376 for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
377 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
378 if (fAlignObjs[iLayer][iModule])
379 delete fAlignObjs[iLayer][iModule];
380 delete [] fAlignObjs[iLayer];
381 }
382 delete [] fAlignObjs;
383}
384
385//______________________________________________________________________________
386void AliAlignmentTracks::Align(Int_t iterations)
387{
388 // This method is just an example
389 // how one can user AlignLayer and
390 // AlignVolume methods to construct
391 // a custom alignment procedure
392 while (iterations > 0) {
393 // First inward pass
394 AlignLayer(AliAlignObj::kTPC1);
395 AlignLayer(AliAlignObj::kSSD2);
396 AlignLayer(AliAlignObj::kSSD1);
397 AlignLayer(AliAlignObj::kSDD2);
398 AlignLayer(AliAlignObj::kSDD1);
399 AlignLayer(AliAlignObj::kSPD2);
400 AlignLayer(AliAlignObj::kSPD1);
401 // Outward pass
402 AlignLayer(AliAlignObj::kSPD2);
403 AlignLayer(AliAlignObj::kSDD1);
404 AlignLayer(AliAlignObj::kSDD2);
405 AlignLayer(AliAlignObj::kSSD1);
406 AlignLayer(AliAlignObj::kSSD2);
407 AlignLayer(AliAlignObj::kTPC1);
408 AlignLayer(AliAlignObj::kTPC2);
409 AlignLayer(AliAlignObj::kTRD1);
410 AlignLayer(AliAlignObj::kTRD2);
411 AlignLayer(AliAlignObj::kTRD3);
412 AlignLayer(AliAlignObj::kTRD4);
413 AlignLayer(AliAlignObj::kTRD5);
414 AlignLayer(AliAlignObj::kTRD6);
415 AlignLayer(AliAlignObj::kTOF);
416 // Again inward
417 AlignLayer(AliAlignObj::kTRD6);
418 AlignLayer(AliAlignObj::kTRD5);
419 AlignLayer(AliAlignObj::kTRD4);
420 AlignLayer(AliAlignObj::kTRD3);
421 AlignLayer(AliAlignObj::kTRD2);
422 AlignLayer(AliAlignObj::kTRD1);
423 AlignLayer(AliAlignObj::kTPC2);
424 }
425}
426
427//______________________________________________________________________________
428void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
429 AliAlignObj::ELayerID layerRangeMin,
430 AliAlignObj::ELayerID layerRangeMax,
431 Int_t iterations)
432{
433 // Align detector volumes within
434 // a given layer.
435 // Tracks are fitted only within
436 // the range defined by the user.
437 while (iterations > 0) {
438 for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
439 UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
440 AlignVolume(volId,layerRangeMin,layerRangeMax);
441 }
442 iterations--;
443 }
444}
445
446//______________________________________________________________________________
447void AliAlignmentTracks::AlignVolume(UShort_t volid,
448 AliAlignObj::ELayerID layerRangeMin,
449 AliAlignObj::ELayerID layerRangeMax)
450{
451 // Align a single detector volume.
452 // Tracks are fitted only within
453 // the range defined by the user.
454
455 // First take the alignment object to be updated
456 Int_t iModule;
457 AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
458 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
459
460 // Then load only the tracks with at least one
461 // space point in the volume (volid)
462 BuildIndex();
463 AliTrackPointArray **points;
464 Int_t nArrays = LoadPoints(volid, points);
465 if (nArrays == 0) return;
466
467 AliTrackResiduals *minimizer = CreateMinimizer();
468 minimizer->SetNTracks(nArrays);
469 minimizer->SetAlignObj(alignObj);
470 AliTrackFitter *fitter = CreateFitter();
471 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
472 fitter->SetTrackPointArray(points[iArray], kFALSE);
473 AliTrackPointArray *pVolId = 0, *pTrack = 0;
474 fitter->Fit(volid,pVolId,pTrack,layerRangeMin,layerRangeMax);
475 minimizer->AddTrackPointArrays(pVolId,pTrack);
476 }
477 minimizer->Minimize();
478
479 UnloadPoints(nArrays, points);
480}
481
482//______________________________________________________________________________
483Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &points)
484{
485 // Load track point arrays with at least
486 // one space point in a given detector
487 // volume (volid).
488 // Use the already created tree index for
489 // fast access.
490 Int_t iModule;
491 AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
492 Int_t nArrays = fLastIndex[iLayer][iModule];
493
494 // In case of empty index
495 if (nArrays == 0) {
496 points = 0;
497 return 0;
498 }
499
500 if (!fPointsTree) {
501 AliWarning("Tree with the space point arrays not initialized!");
502 points = 0;
503 return 0;
504 }
505
506 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
507 TGeoHMatrix m;
508 alignObj->GetMatrix(m);
509 Double_t *rot = m.GetRotationMatrix();
510 Double_t *tr = m.GetTranslation();
511
512 AliTrackPointArray* array = 0;
513 fPointsTree->SetBranchAddress("SP", &array);
514
515 points = new AliTrackPointArray*[nArrays];
516 TArrayI *index = fArrayIndex[iLayer][iModule];
517 AliTrackPoint p;
518 Float_t xyz[3],cov[6];
519 Float_t newxyz[3];
520 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
521 fPointsTree->GetEvent((*index)[iArray]);
522 if (!array) {
523 AliWarning("Wrong space point array index!");
524 continue;
525 }
526 Int_t nPoints = array->GetNPoints();
527 points[iArray] = new AliTrackPointArray(nPoints);
528 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
529 array->GetPoint(p,iPoint);
530 p.GetXYZ(xyz,cov);
531 for (Int_t i = 0; i < 3; i++)
532 newxyz[i] = tr[i]
533 + xyz[0]*rot[3*i]
534 + xyz[1]*rot[3*i+1]
535 + xyz[2]*rot[3*i+2];
536 p.SetXYZ(newxyz,cov);
537 points[iArray]->AddPoint(iPoint,&p);
538 }
539 }
540
541 return nArrays;
542}
543
544//______________________________________________________________________________
545void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
546{
547 // Unload track point arrays for a given
548 // detector volume
549 for (Int_t iArray = 0; iArray < n; iArray++)
550 delete points[iArray];
551 delete [] points;
552}
553
554//______________________________________________________________________________
555AliTrackFitter *AliAlignmentTracks::CreateFitter()
556{
557 // Check if the user has already supplied
558 // a track fitter object.
559 // If not, create a default one.
560 if (!fTrackFitter)
561 fTrackFitter = new AliTrackFitterRieman;
562
563 return fTrackFitter;
564}
565
566//______________________________________________________________________________
567AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
568{
569 // Check if the user has already supplied
570 // a track residuals minimizer object.
571 // If not, create a default one.
572 if (!fMinimizer)
573 fMinimizer = new AliTrackResidualsChi2;
574
575 return fMinimizer;
576}