]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/AliAlignmentTracks.cxx
Switch off branches by hand.
[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 "AliAlignObjParams.h"
31 #include "AliTrackFitterRieman.h"
32 #include "AliTrackResidualsChi2.h"
33 #include "AliESDEvent.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   fAlignObjs(0),
48   fMisalignObjs(0),
49   fTrackFitter(0),
50   fMinimizer(0),
51   fDoUpdate(kTRUE)
52 {
53   // Default constructor
54   InitIndex();
55   InitAlignObjs();
56 }
57
58 //______________________________________________________________________________
59 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
60   fESDChain(esdchain),
61   fPointsFilename("AliTrackPoints.root"),
62   fPointsFile(0),
63   fPointsTree(0),
64   fLastIndex(0),
65   fArrayIndex(0),
66   fIsIndexBuilt(kFALSE),
67   fAlignObjs(0),
68   fMisalignObjs(0),
69   fTrackFitter(0),
70   fMinimizer(0),
71   fDoUpdate(kTRUE)
72 {
73   // Constructor in the case
74   // the user provides an already
75   // built TChain with ESD trees
76   InitIndex();
77   InitAlignObjs();
78 }
79
80
81 //______________________________________________________________________________
82 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
83   fESDChain(new TChain(esdtreename)),
84   fPointsFilename("AliTrackPoints.root"),
85   fPointsFile(0),
86   fPointsTree(0),
87   fLastIndex(0),
88   fArrayIndex(0),
89   fIsIndexBuilt(kFALSE),
90   fAlignObjs(0),
91   fMisalignObjs(0),
92   fTrackFitter(0),
93   fMinimizer(0),
94   fDoUpdate(kTRUE)
95 {
96   // Constructor in the case
97   // the user provides a single ESD file
98   // or a directory containing ESD files
99   fESDChain->Add(esdfilename);
100
101   InitIndex();
102   InitAlignObjs();
103 }
104
105
106 //______________________________________________________________________________
107 AliAlignmentTracks::~AliAlignmentTracks()
108 {
109   // Destructor
110   if (fESDChain) delete fESDChain;
111
112   DeleteIndex();
113   DeleteAlignObjs();
114
115   delete fTrackFitter;
116   delete fMinimizer;
117
118   if (fPointsFile) fPointsFile->Close();
119 }
120
121 //______________________________________________________________________________
122 void AliAlignmentTracks::AddESD(TChain *esdchain)
123 {
124   // Add a chain with ESD files
125   if (fESDChain)
126     fESDChain->Add(esdchain);
127   else
128     fESDChain = esdchain;
129 }
130
131 //______________________________________________________________________________
132 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
133 {
134   // Add a single file or
135   // a directory to the chain
136   // with the ESD files
137   if (fESDChain)
138     fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
139   else {
140     fESDChain = new TChain(esdtreename);
141     fESDChain->Add(esdfilename);
142   }
143 }
144
145 //________________________________________________________________________
146
147 void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
148                                     Int_t minITSpts,
149                                     Bool_t cuts,
150                                     Float_t minMom,Float_t maxMom,
151                                     Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
152                                     Float_t minSinTheta,Float_t maxSinTheta)
153 {
154   // Analyzes and filters ESD tracks
155   // Stores the selected track space points
156   // into the output file
157
158   if (!fESDChain) return;
159
160   AliESDEvent *esd = new AliESDEvent();
161   esd->ReadFromTree(fESDChain);
162   AliESDfriend *esdf = 0; 
163   fESDChain->SetBranchStatus("ESDfriend*",1);
164   fESDChain->SetBranchAddress("ESDfriend.",&esdf);
165
166   // Open the output file
167   if (fPointsFilename.Data() == "") {
168     AliWarning("Incorrect output filename!");
169     return;
170   }
171
172   TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
173   if (!pointsFile || !pointsFile->IsOpen()) {
174     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
175     return;
176   }
177
178   TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
179   const AliTrackPointArray *array = 0;
180   AliTrackPointArray *array2 = 0;
181   if(onlyITS) {   // only ITS AliTrackPoints 
182     pointsTree->Branch("SP","AliTrackPointArray", &array2);
183   } else {
184     pointsTree->Branch("SP","AliTrackPointArray", &array);
185   } 
186
187   Int_t ievent = 0;
188   while (fESDChain->GetEntry(ievent++)) {
189     if (!esd) break;
190
191     esd->SetESDfriend(esdf); //Attach the friend to the ESD
192
193     Int_t ntracks = esd->GetNumberOfTracks();
194     for (Int_t itrack=0; itrack < ntracks; itrack++) {
195       AliESDtrack * track = esd->GetTrack(itrack);
196       if (!track) continue;
197
198       if(track->GetNcls(0) < minITSpts) continue;
199       if(cuts) {
200         if(track->GetP()<minMom || track->GetP()>maxMom) continue;
201         Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
202         if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
203         Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
204         if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
205       } 
206       // UInt_t status = AliESDtrack::kITSpid; 
207       // status|=AliESDtrack::kTPCpid; 
208       // status|=AliESDtrack::kTRDpid; 
209       // if ((track->GetStatus() & status) != status) continue;
210
211       AliTrackPoint point;
212       array = track->GetTrackPointArray();
213
214       if(onlyITS) {
215         array2 = new AliTrackPointArray(track->GetNcls(0));
216         Int_t ipt,volId,modId,layerId;
217         Int_t jpt=0;
218         for(ipt=0; ipt<array->GetNPoints(); ipt++) {
219           array->GetPoint(point,ipt);
220           volId = point.GetVolumeID();
221           layerId = AliGeomManager::VolUIDToLayer(volId,modId);
222           if(layerId<7) {
223             array2->AddPoint(jpt,&point);
224             jpt++;
225           }
226         }
227       } // end if(onlyITS)
228  
229       pointsTree->Fill();
230     }
231   }
232
233   if (!pointsTree->Write()) {
234     AliWarning("Can't write the tree with track point arrays!");
235     return;
236   }
237
238   pointsFile->Close();
239 }
240
241 //______________________________________________________________________________
242 void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
243                                    Int_t minITSpts,Float_t maxMatchingAngle,
244                                    Bool_t cuts,
245                                    Float_t minMom,Float_t maxMom,
246                                    Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
247                                    Float_t minSinTheta,Float_t maxSinTheta)
248 {
249   // Analyzes and filters ESD tracks
250   // Merges inward and outward tracks in one single track
251   // Stores the selected track space points
252   // into the output file
253
254   if (!fESDChain) return;
255
256   AliESDEvent *esd = new AliESDEvent();
257   esd->ReadFromTree(fESDChain);
258   AliESDfriend *esdf = 0; 
259   fESDChain->SetBranchStatus("ESDfriend*",1);
260   fESDChain->SetBranchAddress("ESDfriend.",&esdf);
261
262   // Open the output file
263   if (fPointsFilename.Data() == "") {
264     AliWarning("Incorrect output filename!");
265     return;
266   }
267
268   TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
269   if (!pointsFile || !pointsFile->IsOpen()) {
270     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
271     return;
272   }
273
274   TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
275   const AliTrackPointArray *array = 0;
276   AliTrackPointArray *array2 = 0;
277   pointsTree->Branch("SP","AliTrackPointArray", &array2);
278
279   Int_t ievent = 0;
280   while (fESDChain->GetEntry(ievent++)) {
281     if (!esd) break;
282
283     esd->SetESDfriend(esdf); //Attach the friend to the ESD
284
285     Int_t ntracks = esd->GetNumberOfTracks();
286     if(ntracks<2) continue;
287     Int_t *goodtracksArray = new Int_t[ntracks];
288     Float_t *phiArray = new Float_t[ntracks];
289     Float_t *thetaArray = new Float_t[ntracks];
290     Int_t ngt=0;
291     for (Int_t itrack=0; itrack < ntracks; itrack++) {
292       AliESDtrack * track = esd->GetTrack(itrack);
293       if (!track) continue;
294
295       if(track->GetNcls(0) < minITSpts) continue;
296       Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
297       Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
298       if(cuts) {
299         if(track->GetP()<minMom || track->GetP()>maxMom) continue;
300         Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
301         if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
302         Float_t sintheta = TMath::Sin(theta);
303         if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
304       } 
305       goodtracksArray[ngt]=itrack;
306       phiArray[ngt]=phi;
307       thetaArray[ngt]=theta;
308       ngt++;
309     }
310
311     if(ngt<2) {
312       delete [] goodtracksArray; goodtracksArray=0;
313       delete [] phiArray; phiArray=0;
314       delete [] thetaArray; thetaArray=0;
315       continue;
316     }
317
318     // check matching of the two tracks from the muon
319     Float_t min = 10000000.;
320     Int_t good1 = -1, good2 = -1;
321     for(Int_t itr1=0; itr1<ngt-1; itr1++) {
322       for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
323         Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
324         if(deltatheta>maxMatchingAngle) continue;
325         Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
326         if(deltaphi>maxMatchingAngle) continue;
327         printf("%f  %f     %f  %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
328         if(deltatheta+deltaphi<min) {
329           min=deltatheta+deltaphi;
330           good1 = goodtracksArray[itr1];
331           good2 = goodtracksArray[itr2];
332         }
333       }
334     }
335
336     delete [] goodtracksArray; goodtracksArray=0;
337     delete [] phiArray; phiArray=0;
338     delete [] thetaArray; thetaArray=0;
339
340     if(good1<0) continue;
341
342     AliESDtrack * track1 = esd->GetTrack(good1);
343     AliESDtrack * track2 = esd->GetTrack(good2);
344
345     Int_t ntotpts;
346     if(onlyITS) {
347       ntotpts = track1->GetNcls(0)+track2->GetNcls(0);
348     } else {
349       ntotpts = track1->GetTrackPointArray()->GetNPoints()+
350                 track2->GetTrackPointArray()->GetNPoints();
351     }
352     array2 = new AliTrackPointArray(ntotpts);
353     AliTrackPoint point;
354     Int_t ipt,volId,modId,layerId;
355     Int_t jpt=0;
356     array = track1->GetTrackPointArray();
357     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
358       array->GetPoint(point,ipt);
359       if(onlyITS) {
360         volId = point.GetVolumeID();
361         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
362         if(layerId>6) continue;
363       }
364       array2->AddPoint(jpt,&point);
365       jpt++;
366     }
367     array = track2->GetTrackPointArray();
368     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
369       array->GetPoint(point,ipt);
370       if(onlyITS) {
371         volId = point.GetVolumeID();
372         layerId = AliGeomManager::VolUIDToLayer(volId,modId);
373         if(layerId>6) continue;
374       }
375       array2->AddPoint(jpt,&point);
376       jpt++;
377     }
378
379     pointsTree->Fill();
380   }
381
382   if (!pointsTree->Write()) {
383     AliWarning("Can't write the tree with track point arrays!");
384     return;
385   }
386
387   pointsFile->Close();
388 }
389
390 //______________________________________________________________________________
391 void AliAlignmentTracks::ProcessESD(TSelector *selector)
392 {
393   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
394 }
395
396 //______________________________________________________________________________
397 void AliAlignmentTracks::BuildIndex()
398 {
399   // Build index of points tree entries
400   // Used for access based on the volume IDs
401   if (fIsIndexBuilt) return;
402
403   fIsIndexBuilt = kTRUE;
404
405   // Dummy object is created in order
406   // to initialize the volume paths
407   AliAlignObjParams alobj;
408
409   fPointsFile = TFile::Open(fPointsFilename);
410   if (!fPointsFile || !fPointsFile->IsOpen()) {
411     AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
412     return;
413   }
414   
415   //  AliTrackPointArray* array = new AliTrackPointArray;
416   AliTrackPointArray* array = 0;
417   fPointsTree = (TTree*) fPointsFile->Get("spTree");
418   if (!fPointsTree) {
419     AliWarning("No pointsTree found!");
420     return;
421   }
422   fPointsTree->SetBranchAddress("SP", &array);
423
424   Int_t nArrays = (Int_t)fPointsTree->GetEntries();
425   for (Int_t iArray = 0; iArray < nArrays; iArray++)
426     {
427       fPointsTree->GetEvent(iArray);
428       if (!array) continue;
429       for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
430         UShort_t volId = array->GetVolumeID()[ipoint];
431         // check if the volId is valid
432         if (!AliGeomManager::SymName(volId)) {
433           AliError(Form("The volume id %d has no default volume name !",
434                         volId));
435           continue;
436         }
437         Int_t modId;
438         Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
439                       - AliGeomManager::kFirstLayer;
440         if (!fArrayIndex[layerId][modId]) {
441           //first entry for this volume
442           fArrayIndex[layerId][modId] = new TArrayI(1000);
443         }
444         else {
445           Int_t size = fArrayIndex[layerId][modId]->GetSize();
446           // If needed allocate new size
447           if (fLastIndex[layerId][modId] >= size)
448             fArrayIndex[layerId][modId]->Set(size + 1000);
449         }
450
451         // Check if the index is already filled
452         Bool_t fillIndex = kTRUE;
453         if (fLastIndex[layerId][modId] != 0) {
454           if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
455             fillIndex = kFALSE;
456         }
457         // Fill the index array and store last filled index
458         if (fillIndex) {
459           (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
460           fLastIndex[layerId][modId]++;
461         }
462       }
463     }
464
465 }
466
467 //______________________________________________________________________________
468 void AliAlignmentTracks::InitIndex()
469 {
470   // Initialize the index arrays
471   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
472   fLastIndex = new Int_t*[nLayers];
473   fArrayIndex = new TArrayI**[nLayers];
474   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
475     fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
476     fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
477     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
478       fLastIndex[iLayer][iModule] = 0;
479       fArrayIndex[iLayer][iModule] = 0;
480     }
481   }
482 }
483
484 //______________________________________________________________________________
485 void AliAlignmentTracks::ResetIndex()
486 {
487   // Reset the value of the last filled index
488   // Do not realocate memory
489
490   fIsIndexBuilt = kFALSE;
491   
492   for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
493     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
494       fLastIndex[iLayer][iModule] = 0;
495     }
496   }
497 }
498
499 //______________________________________________________________________________
500 void AliAlignmentTracks::DeleteIndex()
501 {
502   // Delete the index arrays
503   // Called by the destructor
504   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
505     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
506       if (fArrayIndex[iLayer][iModule]) {
507         delete fArrayIndex[iLayer][iModule];
508         fArrayIndex[iLayer][iModule] = 0;
509       }
510     }
511     delete [] fLastIndex[iLayer];
512     delete [] fArrayIndex[iLayer];
513   }
514   delete [] fLastIndex;
515   delete [] fArrayIndex;
516 }
517
518 //______________________________________________________________________________
519 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
520 {
521   // Read alignment object from a file
522   // To be replaced by a call to CDB
523   AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
524
525   return kFALSE;
526 }
527
528 //______________________________________________________________________________
529 void AliAlignmentTracks::InitAlignObjs()
530 {
531   // Initialize the alignment objects array
532   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
533   fAlignObjs = new AliAlignObj**[nLayers];
534   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
535     fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
536     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
537       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
538       fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
539     }
540   }
541 }
542
543 //______________________________________________________________________________
544 void AliAlignmentTracks::ResetAlignObjs()
545 {
546   // Reset the alignment objects array
547   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
548     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
549       fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
550   }
551 }
552
553 //______________________________________________________________________________
554 void AliAlignmentTracks::DeleteAlignObjs()
555 {
556   // Delete the alignment objects array
557   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
558     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
559       if (fAlignObjs[iLayer][iModule])
560         delete fAlignObjs[iLayer][iModule];
561     delete [] fAlignObjs[iLayer];
562   }
563   delete [] fAlignObjs;
564   fAlignObjs = 0;
565 }
566
567 void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
568                                        AliGeomManager::ELayerID lastLayer,
569                                        AliGeomManager::ELayerID layerRangeMin,
570                                        AliGeomManager::ELayerID layerRangeMax,
571                                        Int_t iterations)
572 {
573   // Align detector volumes within
574   // a given layer range
575   // (could be whole detector).
576   // Tracks are fitted only within
577   // the range defined by the user.
578   Int_t nModules = 0;
579   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
580     nModules += AliGeomManager::LayerSize(iLayer);
581   TArrayI volIds(nModules);
582
583   Int_t modnum = 0;
584   for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
585     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
586       UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
587       volIds.AddAt(volId,modnum);
588       modnum++;
589     }
590   }
591
592   while (iterations > 0) {
593     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
594     iterations--;
595   }
596 }
597
598 //______________________________________________________________________________
599 void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
600                                     AliGeomManager::ELayerID layerRangeMin,
601                                     AliGeomManager::ELayerID layerRangeMax,
602                                     Int_t iterations)
603 {
604   // Align detector volumes within
605   // a given layer.
606   // Tracks are fitted only within
607   // the range defined by the user.
608   Int_t nModules = AliGeomManager::LayerSize(layer);
609   TArrayI volIds(nModules);
610   for (Int_t iModule = 0; iModule < nModules; iModule++) {
611     UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
612     volIds.AddAt(volId,iModule);
613   }
614
615   while (iterations > 0) {
616     AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
617     iterations--;
618   }
619 }
620
621 //______________________________________________________________________________
622 void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
623                                      Int_t iterations)
624 {
625   // Align single detector volume to
626   // another volume.
627   // Tracks are fitted only within
628   // the second volume.
629   TArrayI volIds(1);
630   volIds.AddAt(volId,0);
631   TArrayI volIdsFit(1);
632   volIdsFit.AddAt(volIdFit,0);
633
634   while (iterations > 0) {
635     AlignVolumes(&volIds,&volIdsFit);
636     iterations--;
637   }
638 }
639
640 //______________________________________________________________________________
641 void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
642                                      AliGeomManager::ELayerID layerRangeMin,
643                                      AliGeomManager::ELayerID layerRangeMax,
644                                      Int_t iterations)
645 {
646   // Align a set of detector volumes.
647   // Tracks are fitted only within
648   // the range defined by the user
649   // (by layerRangeMin and layerRangeMax)
650   // or within the set of volidsfit
651   // Repeat the procedure 'iterations' times
652
653   Int_t nVolIds = volids->GetSize();
654   if (nVolIds == 0) {
655     AliError("Volume IDs array is empty!");
656     return;
657   }
658
659   // Load only the tracks with at least one
660   // space point in the set of volume (volids)
661   BuildIndex();
662   AliTrackPointArray **points;
663   // Start the iterations
664   while (iterations > 0) {
665     Int_t nArrays = LoadPoints(volids, points);
666     if (nArrays == 0) return;
667
668     AliTrackResiduals *minimizer = CreateMinimizer();
669     minimizer->SetNTracks(nArrays);
670     minimizer->InitAlignObj();
671     AliTrackFitter *fitter = CreateFitter();
672     for (Int_t iArray = 0; iArray < nArrays; iArray++) {
673       if (!points[iArray]) continue;
674       fitter->SetTrackPointArray(points[iArray], kFALSE);
675       if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
676       AliTrackPointArray *pVolId,*pTrack;
677       fitter->GetTrackResiduals(pVolId,pTrack);
678       minimizer->AddTrackPointArrays(pVolId,pTrack);
679     }
680     minimizer->Minimize();
681
682     // Update the alignment object(s)
683     if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
684       UShort_t volid = (*volids)[iVolId];
685       Int_t iModule;
686       AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
687       AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];      
688       *alignObj *= *minimizer->GetAlignObj();
689       alignObj->Print("");
690     }
691
692     UnloadPoints(nArrays, points);
693     
694     iterations--;
695   }
696 }
697   
698 //______________________________________________________________________________
699 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
700 {
701   // Load track point arrays with at least
702   // one space point in a given set of detector
703   // volumes (array volids).
704   // Use the already created tree index for
705   // fast access.
706
707   if (!fPointsTree) {
708     AliError("Tree with the space point arrays not initialized!");
709     points = 0;
710     return 0;
711   }
712
713   Int_t nVolIds = volids->GetSize();
714   if (nVolIds == 0) {
715     AliError("Volume IDs array is empty!");
716     points = 0;
717     return 0;
718   }
719
720   Int_t nArrays = 0;
721   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
722     UShort_t volid = (*volids)[iVolId];
723     Int_t iModule;
724     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
725
726     // In case of empty index
727     if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
728       AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
729       continue;
730     }
731     nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
732   }
733
734   if (nArrays == 0) {
735     AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
736     points = 0;
737     return 0;
738   }
739
740   AliTrackPointArray* array = 0;
741   fPointsTree->SetBranchAddress("SP", &array);
742
743   // Allocate the pointer to the space-point arrays
744   points = new AliTrackPointArray*[nArrays];
745   for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
746
747   // Init the array used to flag already loaded tree entries
748   Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
749   for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
750     indexUsed[i] = kFALSE;
751
752   // Start the loop over the volume ids
753   Int_t iArray = 0;
754   for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
755     UShort_t volid = (*volids)[iVolId];
756     Int_t iModule;
757     AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
758
759     Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
760     TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
761     AliTrackPoint p;
762
763     for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
764
765       // Get tree entry
766       Int_t entry = (*index)[iArrayId];
767       if (indexUsed[entry] == kTRUE) continue;
768       fPointsTree->GetEvent(entry);
769       if (!array) {
770         AliWarning("Wrong space point array index!");
771         continue;
772       }
773       indexUsed[entry] = kTRUE;
774
775       // Get the space-point array
776       Int_t nPoints = array->GetNPoints();
777       points[iArray] = new AliTrackPointArray(nPoints);
778       for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
779         array->GetPoint(p,iPoint);
780         Int_t modnum;
781         AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
782         // check if the layer id is valid
783         if ((layer < AliGeomManager::kFirstLayer) ||
784             (layer >= AliGeomManager::kLastLayer)) {
785           AliError(Form("Layer index is invalid: %d (%d -> %d) !",
786                         layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
787           continue;
788         }
789         if ((modnum >= AliGeomManager::LayerSize(layer)) ||
790             (modnum < 0)) {
791           AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
792                         layer,modnum,AliGeomManager::LayerSize(layer)));
793           continue;
794         }
795
796         // Misalignment is introduced here
797         // Switch it off in case of real
798         // alignment job!
799         if (fMisalignObjs) {
800           AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
801           if (misalignObj)
802             misalignObj->Transform(p);
803         }
804         // End of misalignment
805
806         AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
807         alignObj->Transform(p);
808         points[iArray]->AddPoint(iPoint,&p);
809       }
810       iArray++;
811     }
812   }
813
814
815   delete [] indexUsed;
816
817   return nArrays;
818 }
819
820 //______________________________________________________________________________
821 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
822 {
823   // Unload track point arrays for a given
824   // detector volume
825   for (Int_t iArray = 0; iArray < n; iArray++)
826     delete points[iArray];
827   delete [] points;
828 }
829
830 //______________________________________________________________________________
831 AliTrackFitter *AliAlignmentTracks::CreateFitter()
832 {
833   // Check if the user has already supplied
834   // a track fitter object.
835   // If not, create a default one.
836   if (!fTrackFitter)
837     fTrackFitter = new AliTrackFitterRieman;
838
839   return fTrackFitter;
840 }
841
842 //______________________________________________________________________________
843 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
844 {
845   // Check if the user has already supplied
846   // a track residuals minimizer object.
847   // If not, create a default one.
848   if (!fMinimizer)
849     fMinimizer = new AliTrackResidualsChi2;
850
851   return fMinimizer;
852 }
853
854 //______________________________________________________________________________
855 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
856 {
857   // The method reads from a file a set of AliAlignObj which are
858   // then used to apply misalignments directly on the track
859   // space-points. The method is supposed to be used only for
860   // fast development and debugging of the alignment algorithms.
861   // Be careful not to use it in the case of 'real' alignment
862   // scenario since it will bias the results.
863
864   // Initialize the misalignment objects array
865   Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
866   fMisalignObjs = new AliAlignObj**[nLayers];
867   for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
868     fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
869     for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
870       fMisalignObjs[iLayer][iModule] = 0x0;
871   }
872
873   // Open the misliagnment file and load the array with
874   // misalignment objects
875   TFile* inFile = TFile::Open(misalignObjFileName,"READ");
876   if (!inFile || !inFile->IsOpen()) {
877     AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
878     return kFALSE;
879   }
880
881   TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
882   if (!array) {
883     AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
884     inFile->Close();
885     return kFALSE;
886   }
887   inFile->Close();
888
889   // Store the misalignment objects for further usage  
890   Int_t nObjs = array->GetEntriesFast();
891   AliGeomManager::ELayerID layerId; // volume layer
892   Int_t modId; // volume ID inside the layer
893   for(Int_t i=0; i<nObjs; i++)
894     {
895       AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
896       alObj->GetVolUID(layerId,modId);
897       fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
898     }
899   return kTRUE;
900 }