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