AliRoot Core  3dc7879 (3dc7879)
AliAlignmentTracks.cxx
Go to the documentation of this file.
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //_____________________________________________________________________________
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 //______________________________________________________________________________
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 
486  if (!fPointsFile || !fPointsFile->IsOpen()) {
487  AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
488  return;
489  }
490 
491  // AliTrackPointArray* array = new AliTrackPointArray;
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)
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 //______________________________________________________________________________
545 {
546  // Initialize the index arrays
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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;
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 //______________________________________________________________________________
632 {
633  // Initialize the alignment objects array
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 
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 //______________________________________________________________________________
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;
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;
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 
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;
876 
877  Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
878  TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
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;
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) !",
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 //______________________________________________________________________________
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 //______________________________________________________________________________
966 {
967  // Check if the user has already supplied
968  // a track fitter object.
969  // If not, create a default one.
970  if (!fTrackFitter)
972 
973  return fTrackFitter;
974 }
975 
976 //______________________________________________________________________________
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)
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
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 //________________________________________________
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 }
Bool_t AddPoint(Int_t i, const AliTrackPoint *p)
Int_t GetNcls(Int_t idet) const
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TFile * Open(const char *filename, Long64_t &nevents)
AliTrackResiduals * CreateMinimizer()
Int_t LoadPoints(const TArrayI *volids, AliTrackPointArray **&points, Int_t &pointsdim)
void AddESD(TChain *esdchain)
TString fPointsFilename
Chain with ESDs.
static const char * SymName(UShort_t voluid)
UShort_t GetVolumeID() const
void Print(Option_t *) const
Bool_t Misalign(const char *misalignObjFileName, const char *arrayName)
Int_t GetNPoints() const
void Transform(AliTrackPoint &p, Bool_t copycov=kFALSE) const
AliAlignObj *** fMisalignObjs
static Bool_t GetOrigRotation(Int_t index, Double_t r[9])
Float_t p[]
Definition: kNNTest.C:133
AliTrackFitter * CreateFitter()
Double_t GetAlpha() const
void GetTrackResiduals(AliTrackPointArray *&pVolId, AliTrackPointArray *&pTrack) const
static Int_t LayerSize(Int_t layerId)
AliTPCfastTrack * track
virtual Double_t GetTgl() const
#define AliWarning(message)
Definition: AliLog.h:541
const UShort_t * GetVolumeID() const
TObjArray * array
Definition: AnalyzeLaser.C:12
void ProcessESDCosmics(Bool_t onlyITS=kFALSE, Int_t minITSpts=0, Float_t maxMatchingAngle=0.17, Bool_t cuts=kTRUE, Float_t minAngleWrtITSModulePlanes=0., Float_t minMom=0.3, Float_t maxMom=1.e9, Float_t minAbsSinPhi=0., Float_t maxAbsSinPhi=1., Float_t minSinTheta=0., Float_t maxSinTheta=1.)
AliTrackResiduals * fMinimizer
Bool_t AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit=0x0, AliGeomManager::ELayerID layerRangeMin=AliGeomManager::kFirstLayer, AliGeomManager::ELayerID layerRangeMax=AliGeomManager::kLastLayer, Int_t iterations=1)
Bool_t GetPoint(AliTrackPoint &p, Int_t i) const
static ELayerID VolUIDToLayer(UShort_t voluid, Int_t &modId)
const AliTrackPointArray * GetTrackPointArray() const
Definition: AliESDtrack.h:429
Bool_t AlignLayer(AliGeomManager::ELayerID layer, AliGeomManager::ELayerID layerRangeMin=AliGeomManager::kFirstLayer, AliGeomManager::ELayerID layerRangeMax=AliGeomManager::kLastLayer, Int_t iterations=1)
Bool_t AlignDetector(AliGeomManager::ELayerID firstLayer, AliGeomManager::ELayerID lastLayer, AliGeomManager::ELayerID layerRangeMin=AliGeomManager::kFirstLayer, AliGeomManager::ELayerID layerRangeMax=AliGeomManager::kLastLayer, Int_t iterations=1)
void SetESDfriend(const AliESDfriend *f) const
void UnloadPoints(Int_t n, AliTrackPointArray **points)
void ReadFromTree(TTree *tree, Option_t *opt="")
AliAlignObj * GetAlignObj() const
Bool_t fIsIndexBuilt
Volume arrays which contains the tree index.
Bool_t ReadAlignObjs(const char *alignObjFileName="AlignObjs.root", const char *arrayName="Alignment")
virtual Bool_t Minimize()=0
AliAlignObj *** fAlignObjs
Int_t GetNumberOfTracks() const
Definition: AliESDEvent.h:536
AliESDtrack * GetTrack(Int_t i) const
Definition: AliESDEvent.h:405
static UShort_t LayerToVolUID(ELayerID layerId, Int_t modId)
void SetNTracks(Int_t ntracks)
const char * GetSymName() const
Definition: AliAlignObj.h:63
UShort_t GetVolUID() const
Definition: AliAlignObj.h:64
void WriteRealignObjArray(TString outfilename, AliGeomManager::ELayerID layerRangeMin, AliGeomManager::ELayerID layerRangeMax)
#define AliError(message)
Definition: AliLog.h:591
Bool_t AddTrackPointArrays(AliTrackPointArray *volarray, AliTrackPointArray *trackarray)
void ProcessESD(TSelector *selector)
TArrayI *** fArrayIndex
Last filled index in volume arrays.
Bool_t AlignVolume(UShort_t volId, UShort_t volIdFit, Int_t iterations)
void GetDirection(Double_t d[3]) const
virtual void SetTrackPointArray(AliTrackPointArray *array, Bool_t owner=kTRUE)
AliTrackFitter * fTrackFitter
virtual Bool_t Fit(const TArrayI *volIds, const TArrayI *volIdsFit=0x0, AliGeomManager::ELayerID layerRangeMin=AliGeomManager::kFirstLayer, AliGeomManager::ELayerID layerRangeMax=AliGeomManager::kLastLayer)