089decebacc35ccb5a9f38970cab5b0d1bd01e8f
[u/mrichter/AliRoot.git] / STEER / AliAlignmentTracks.cxx
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 #include <TFile.h>
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
36 ClassImp(AliAlignmentTracks)
37
38 //______________________________________________________________________________
39 AliAlignmentTracks::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 //______________________________________________________________________________
56 AliAlignmentTracks::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 //______________________________________________________________________________
76 AliAlignmentTracks::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 //______________________________________________________________________________
97 AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
98   TObject(alignment)
99 {
100   // Copy constructor
101   // not implemented
102   AliWarning("Copy constructor not implemented!");
103 }
104
105 //______________________________________________________________________________
106 AliAlignmentTracks& 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 //______________________________________________________________________________
120 AliAlignmentTracks::~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 //______________________________________________________________________________
135 void 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 //______________________________________________________________________________
145 void 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 //______________________________________________________________________________
159 void 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 //______________________________________________________________________________
215 void AliAlignmentTracks::ProcessESD(TSelector *selector)
216 {
217   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
218 }
219
220 //______________________________________________________________________________
221 void 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 //______________________________________________________________________________
293 void 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 //______________________________________________________________________________
310 void 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 //______________________________________________________________________________
322 void 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 //______________________________________________________________________________
341 Bool_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 //______________________________________________________________________________
351 void 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 //______________________________________________________________________________
364 void 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 //______________________________________________________________________________
374 void 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 //______________________________________________________________________________
387 void 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 //______________________________________________________________________________
429 void 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 //______________________________________________________________________________
448 void 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 //______________________________________________________________________________
484 Int_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 //______________________________________________________________________________
546 void 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 //______________________________________________________________________________
556 AliTrackFitter *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 //______________________________________________________________________________
568 AliTrackResiduals *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 }