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