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