1c8f2b64dfb352b81db8241df6affcf94ee182c0
[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
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
35 ClassImp(AliAlignmentTracks)
36
37 //______________________________________________________________________________
38 AliAlignmentTracks::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 //______________________________________________________________________________
55 AliAlignmentTracks::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 //______________________________________________________________________________
75 AliAlignmentTracks::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 //______________________________________________________________________________
96 AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
97   TObject(alignment)
98 {
99   // Copy constructor
100   // not implemented
101   AliWarning("Copy constructor not implemented!");
102 }
103
104 //______________________________________________________________________________
105 AliAlignmentTracks& 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 //______________________________________________________________________________
119 AliAlignmentTracks::~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 //______________________________________________________________________________
134 void 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 //______________________________________________________________________________
144 void 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 //______________________________________________________________________________
158 void 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 //______________________________________________________________________________
214 void AliAlignmentTracks::ProcessESD(TSelector *selector)
215 {
216   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
217 }
218
219 //______________________________________________________________________________
220 void 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 //______________________________________________________________________________
292 void 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 //______________________________________________________________________________
309 void 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 //______________________________________________________________________________
321 void 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 //______________________________________________________________________________
340 Bool_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 //______________________________________________________________________________
350 void 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 //______________________________________________________________________________
363 void 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 //______________________________________________________________________________
373 void 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 //______________________________________________________________________________
386 void 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 //______________________________________________________________________________
428 void 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 //______________________________________________________________________________
447 void 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 //______________________________________________________________________________
483 Int_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 //______________________________________________________________________________
545 void 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 //______________________________________________________________________________
555 AliTrackFitter *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 //______________________________________________________________________________
567 AliTrackResiduals *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 }