1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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 //-----------------------------------------------------------------
29 #include "AliAlignmentTracks.h"
30 #include "AliTrackPointArray.h"
31 #include "AliAlignObjParams.h"
32 #include "AliTrackFitterRieman.h"
33 #include "AliTrackResidualsChi2.h"
34 #include "AliESDEvent.h"
37 ClassImp(AliAlignmentTracks)
39 //______________________________________________________________________________
40 AliAlignmentTracks::AliAlignmentTracks():
42 fPointsFilename("AliTrackPoints.root"),
47 fIsIndexBuilt(kFALSE),
54 // Default constructor
59 //______________________________________________________________________________
60 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
62 fPointsFilename("AliTrackPoints.root"),
67 fIsIndexBuilt(kFALSE),
74 // Constructor in the case
75 // the user provides an already
76 // built TChain with ESD trees
82 //______________________________________________________________________________
83 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
84 fESDChain(new TChain(esdtreename)),
85 fPointsFilename("AliTrackPoints.root"),
90 fIsIndexBuilt(kFALSE),
97 // Constructor in the case
98 // the user provides a single ESD file
99 // or a directory containing ESD files
100 fESDChain->Add(esdfilename);
107 //______________________________________________________________________________
108 AliAlignmentTracks::~AliAlignmentTracks()
111 if (fESDChain) delete fESDChain;
119 if (fPointsFile) fPointsFile->Close();
122 //______________________________________________________________________________
123 void AliAlignmentTracks::AddESD(TChain *esdchain)
125 // Add a chain with ESD files
127 fESDChain->Add(esdchain);
129 fESDChain = esdchain;
132 //______________________________________________________________________________
133 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
135 // Add a single file or
136 // a directory to the chain
137 // with the ESD files
139 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
141 fESDChain = new TChain(esdtreename);
142 fESDChain->Add(esdfilename);
147 //________________________________________________________________________
148 void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
151 Float_t minAngleWrtITSModulePlanes,
152 Float_t minMom,Float_t maxMom,
153 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
154 Float_t minSinTheta,Float_t maxSinTheta)
156 // Analyzes and filters ESD tracks
157 // Stores the selected track space points
158 // into the output file
160 if (!fESDChain) return;
162 AliESDEvent *esd = new AliESDEvent();
163 esd->ReadFromTree(fESDChain);
164 AliESDfriend *esdf = 0;
165 fESDChain->SetBranchStatus("ESDfriend*",1);
166 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
168 // Open the output file
169 if (fPointsFilename.Data() == "") {
170 AliWarning("Incorrect output filename!");
174 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
175 if (!pointsFile || !pointsFile->IsOpen()) {
176 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
180 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
181 const AliTrackPointArray *array = 0;
182 AliTrackPointArray *array2 = 0;
183 if(onlyITS) { // only ITS AliTrackPoints
184 pointsTree->Branch("SP","AliTrackPointArray", &array2);
186 pointsTree->Branch("SP","AliTrackPointArray", &array);
190 while (fESDChain->GetEntry(ievent++)) {
193 esd->SetESDfriend(esdf); //Attach the friend to the ESD
195 Int_t ntracks = esd->GetNumberOfTracks();
196 for (Int_t itrack=0; itrack < ntracks; itrack++) {
197 AliESDtrack * track = esd->GetTrack(itrack);
198 if (!track) continue;
200 if(track->GetNcls(0) < minITSpts) continue;
202 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
203 Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
204 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
205 Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
206 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
210 array = track->GetTrackPointArray();
213 Bool_t layerOK[6]={kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
214 Int_t ipt,volId,modId,layerId;
216 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
217 array->GetPoint(point,ipt);
218 volId = point.GetVolumeID();
219 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
220 if(layerId>6) continue;
221 // check minAngleWrtITSModulePlanes
223 Double_t p[3]; track->GetDirection(p);
224 TVector3 pvec(p[0],p[1],p[2]);
225 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
226 TVector3 normvec(rot[1],rot[4],rot[7]);
227 Double_t angle = pvec.Angle(normvec);
228 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
229 angle = 0.5*TMath::Pi()-angle;
230 if(angle<minAngleWrtITSModulePlanes) {
231 layerOK[layerId-1]=kFALSE;
237 if(jpt < minITSpts) continue;
238 array2 = new AliTrackPointArray(jpt);
240 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
241 array->GetPoint(point,ipt);
242 volId = point.GetVolumeID();
243 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
244 if(layerId>6 || !layerOK[layerId-1]) continue;
245 array2->AddPoint(jpt,&point);
254 if (!pointsTree->Write()) {
255 AliWarning("Can't write the tree with track point arrays!");
264 //_____________________________________________________________________________
265 void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
266 Int_t minITSpts,Float_t maxMatchingAngle,
268 Float_t minAngleWrtITSModulePlanes,
269 Float_t minMom,Float_t maxMom,
270 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
271 Float_t minSinTheta,Float_t maxSinTheta)
273 // Analyzes and filters ESD tracks
274 // Merges inward and outward tracks in one single track
275 // Stores the selected track space points
276 // into the output file
278 if (!fESDChain) return;
280 AliESDEvent *esd = new AliESDEvent();
281 esd->ReadFromTree(fESDChain);
282 AliESDfriend *esdf = 0;
283 fESDChain->SetBranchStatus("ESDfriend*",1);
284 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
286 // Open the output file
287 if (fPointsFilename.Data() == "") {
288 AliWarning("Incorrect output filename!");
292 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
293 if (!pointsFile || !pointsFile->IsOpen()) {
294 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
298 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
299 const AliTrackPointArray *array = 0;
300 AliTrackPointArray *array2 = 0;
301 pointsTree->Branch("SP","AliTrackPointArray", &array2);
304 while (fESDChain->GetEntry(ievent++)) {
307 esd->SetESDfriend(esdf); //Attach the friend to the ESD
309 Int_t ntracks = esd->GetNumberOfTracks();
310 if(ntracks<2) continue;
311 Int_t *goodtracksArray = new Int_t[ntracks];
312 Float_t *phiArray = new Float_t[ntracks];
313 Float_t *thetaArray = new Float_t[ntracks];
315 for (Int_t itrack=0; itrack < ntracks; itrack++) {
316 AliESDtrack * track = esd->GetTrack(itrack);
317 if (!track) continue;
319 if(track->GetNcls(0) < minITSpts) continue;
320 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
321 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
323 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
324 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
325 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
326 Float_t sintheta = TMath::Sin(theta);
327 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
329 goodtracksArray[ngt]=itrack;
331 thetaArray[ngt]=theta;
336 delete [] goodtracksArray; goodtracksArray=0;
337 delete [] phiArray; phiArray=0;
338 delete [] thetaArray; thetaArray=0;
342 // check matching of the two tracks from the muon
343 Float_t min = 10000000.;
344 Int_t good1 = -1, good2 = -1;
345 for(Int_t itr1=0; itr1<ngt-1; itr1++) {
346 for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
347 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
348 if(deltatheta>maxMatchingAngle) continue;
349 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
350 if(deltaphi>maxMatchingAngle) continue;
351 //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
352 if(deltatheta+deltaphi<min) {
353 min=deltatheta+deltaphi;
354 good1 = goodtracksArray[itr1];
355 good2 = goodtracksArray[itr2];
360 delete [] goodtracksArray; goodtracksArray=0;
361 delete [] phiArray; phiArray=0;
362 delete [] thetaArray; thetaArray=0;
364 if(good1<0) continue;
366 AliESDtrack * track1 = esd->GetTrack(good1);
367 AliESDtrack * track2 = esd->GetTrack(good2);
370 Int_t ipt,volId,modId,layerId;
372 Bool_t layerOK[6][2];
373 for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kTRUE;
374 array = track1->GetTrackPointArray();
375 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
376 array->GetPoint(point,ipt);
378 volId = point.GetVolumeID();
379 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
380 if(layerId>6) continue;
381 // check minAngleWrtITSModulePlanes
383 Double_t p[3]; track1->GetDirection(p);
384 TVector3 pvec(p[0],p[1],p[2]);
385 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
386 TVector3 normvec(rot[1],rot[4],rot[7]);
387 Double_t angle = pvec.Angle(normvec);
388 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
389 angle = 0.5*TMath::Pi()-angle;
390 if(angle<minAngleWrtITSModulePlanes) {
391 layerOK[layerId-1][0]=kFALSE;
398 array = track2->GetTrackPointArray();
399 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
400 array->GetPoint(point,ipt);
402 volId = point.GetVolumeID();
403 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
404 if(layerId>6) continue;
405 // check minAngleWrtITSModulePlanes
407 Double_t p[3]; track2->GetDirection(p);
408 TVector3 pvec(p[0],p[1],p[2]);
409 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
410 TVector3 normvec(rot[1],rot[4],rot[7]);
411 Double_t angle = pvec.Angle(normvec);
412 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
413 angle = 0.5*TMath::Pi()-angle;
414 if(angle<minAngleWrtITSModulePlanes) {
415 layerOK[layerId-1][0]=kFALSE;
423 if(jpt < 2*minITSpts) continue;
424 array2 = new AliTrackPointArray(jpt);
426 array = track1->GetTrackPointArray();
427 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
428 array->GetPoint(point,ipt);
430 volId = point.GetVolumeID();
431 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
432 if(layerId>6 || !layerOK[layerId-1][0]) continue;
434 array2->AddPoint(jpt,&point);
437 array = track2->GetTrackPointArray();
438 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
439 array->GetPoint(point,ipt);
441 volId = point.GetVolumeID();
442 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
443 if(layerId>6 || !layerOK[layerId-1][1]) continue;
445 array2->AddPoint(jpt,&point);
452 if (!pointsTree->Write()) {
453 AliWarning("Can't write the tree with track point arrays!");
461 //______________________________________________________________________________
462 void AliAlignmentTracks::ProcessESD(TSelector *selector)
464 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
467 //______________________________________________________________________________
468 void AliAlignmentTracks::BuildIndex()
470 // Build index of points tree entries
471 // Used for access based on the volume IDs
472 if (fIsIndexBuilt) return;
474 fIsIndexBuilt = kTRUE;
476 // Dummy object is created in order
477 // to initialize the volume paths
478 AliAlignObjParams alobj;
480 fPointsFile = TFile::Open(fPointsFilename);
481 if (!fPointsFile || !fPointsFile->IsOpen()) {
482 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
486 // AliTrackPointArray* array = new AliTrackPointArray;
487 AliTrackPointArray* array = 0;
488 fPointsTree = (TTree*) fPointsFile->Get("spTree");
490 AliWarning("No pointsTree found!");
493 fPointsTree->SetBranchAddress("SP", &array);
495 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
496 for (Int_t iArray = 0; iArray < nArrays; iArray++)
498 fPointsTree->GetEvent(iArray);
499 if (!array) continue;
500 for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
501 UShort_t volId = array->GetVolumeID()[ipoint];
502 // check if the volId is valid
503 if (!AliGeomManager::SymName(volId)) {
504 AliError(Form("The volume id %d has no default volume name !",
509 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
510 - AliGeomManager::kFirstLayer;
511 if (!fArrayIndex[layerId][modId]) {
512 //first entry for this volume
513 fArrayIndex[layerId][modId] = new TArrayI(1000);
516 Int_t size = fArrayIndex[layerId][modId]->GetSize();
517 // If needed allocate new size
518 if (fLastIndex[layerId][modId] >= size)
519 fArrayIndex[layerId][modId]->Set(size + 1000);
522 // Check if the index is already filled
523 Bool_t fillIndex = kTRUE;
524 if (fLastIndex[layerId][modId] != 0) {
525 if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
528 // Fill the index array and store last filled index
530 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
531 fLastIndex[layerId][modId]++;
538 //______________________________________________________________________________
539 void AliAlignmentTracks::InitIndex()
541 // Initialize the index arrays
542 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
543 fLastIndex = new Int_t*[nLayers];
544 fArrayIndex = new TArrayI**[nLayers];
545 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
546 fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
547 fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
548 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
549 fLastIndex[iLayer][iModule] = 0;
550 fArrayIndex[iLayer][iModule] = 0;
555 //______________________________________________________________________________
556 void AliAlignmentTracks::ResetIndex()
558 // Reset the value of the last filled index
559 // Do not realocate memory
561 fIsIndexBuilt = kFALSE;
563 for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
564 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
565 fLastIndex[iLayer][iModule] = 0;
570 //______________________________________________________________________________
571 void AliAlignmentTracks::DeleteIndex()
573 // Delete the index arrays
574 // Called by the destructor
575 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
576 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
577 if (fArrayIndex[iLayer][iModule]) {
578 delete fArrayIndex[iLayer][iModule];
579 fArrayIndex[iLayer][iModule] = 0;
582 delete [] fLastIndex[iLayer];
583 delete [] fArrayIndex[iLayer];
585 delete [] fLastIndex;
586 delete [] fArrayIndex;
589 //______________________________________________________________________________
590 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
592 // Read alignment object from a file
593 // To be replaced by a call to CDB
594 AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
599 //______________________________________________________________________________
600 void AliAlignmentTracks::InitAlignObjs()
602 // Initialize the alignment objects array
603 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
604 fAlignObjs = new AliAlignObj**[nLayers];
605 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
606 fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
607 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
608 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
609 fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
614 //______________________________________________________________________________
615 void AliAlignmentTracks::ResetAlignObjs()
617 // Reset the alignment objects array
618 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
619 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
620 fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
624 //______________________________________________________________________________
625 void AliAlignmentTracks::DeleteAlignObjs()
627 // Delete the alignment objects array
628 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
629 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
630 if (fAlignObjs[iLayer][iModule])
631 delete fAlignObjs[iLayer][iModule];
632 delete [] fAlignObjs[iLayer];
634 delete [] fAlignObjs;
638 void AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
639 AliGeomManager::ELayerID lastLayer,
640 AliGeomManager::ELayerID layerRangeMin,
641 AliGeomManager::ELayerID layerRangeMax,
644 // Align detector volumes within
645 // a given layer range
646 // (could be whole detector).
647 // Tracks are fitted only within
648 // the range defined by the user.
650 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
651 nModules += AliGeomManager::LayerSize(iLayer);
652 TArrayI volIds(nModules);
655 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
656 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
657 UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
658 volIds.AddAt(volId,modnum);
663 while (iterations > 0) {
664 AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
669 //______________________________________________________________________________
670 void AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
671 AliGeomManager::ELayerID layerRangeMin,
672 AliGeomManager::ELayerID layerRangeMax,
675 // Align detector volumes within
677 // Tracks are fitted only within
678 // the range defined by the user.
679 Int_t nModules = AliGeomManager::LayerSize(layer);
680 TArrayI volIds(nModules);
681 for (Int_t iModule = 0; iModule < nModules; iModule++) {
682 UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
683 volIds.AddAt(volId,iModule);
686 while (iterations > 0) {
687 AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax);
692 //______________________________________________________________________________
693 void AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
696 // Align single detector volume to
698 // Tracks are fitted only within
699 // the second volume.
701 volIds.AddAt(volId,0);
702 TArrayI volIdsFit(1);
703 volIdsFit.AddAt(volIdFit,0);
705 while (iterations > 0) {
706 AlignVolumes(&volIds,&volIdsFit);
711 //______________________________________________________________________________
712 void AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
713 AliGeomManager::ELayerID layerRangeMin,
714 AliGeomManager::ELayerID layerRangeMax,
717 // Align a set of detector volumes.
718 // Tracks are fitted only within
719 // the range defined by the user
720 // (by layerRangeMin and layerRangeMax)
721 // or within the set of volidsfit
722 // Repeat the procedure 'iterations' times
724 Int_t nVolIds = volids->GetSize();
726 AliError("Volume IDs array is empty!");
730 // Load only the tracks with at least one
731 // space point in the set of volume (volids)
733 AliTrackPointArray **points;
734 // Start the iterations
735 while (iterations > 0) {
736 Int_t nArrays = LoadPoints(volids, points);
737 if (nArrays == 0) return;
739 AliTrackResiduals *minimizer = CreateMinimizer();
740 minimizer->SetNTracks(nArrays);
741 minimizer->InitAlignObj();
742 AliTrackFitter *fitter = CreateFitter();
743 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
744 if (!points[iArray]) continue;
745 fitter->SetTrackPointArray(points[iArray], kFALSE);
746 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
747 AliTrackPointArray *pVolId,*pTrack;
748 fitter->GetTrackResiduals(pVolId,pTrack);
749 minimizer->AddTrackPointArrays(pVolId,pTrack);
751 minimizer->Minimize();
753 // Update the alignment object(s)
754 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
755 UShort_t volid = (*volids)[iVolId];
757 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
758 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
759 *alignObj *= *minimizer->GetAlignObj();
760 if(iterations==1)alignObj->Print("");
763 UnloadPoints(nArrays, points);
769 //______________________________________________________________________________
770 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points)
772 // Load track point arrays with at least
773 // one space point in a given set of detector
774 // volumes (array volids).
775 // Use the already created tree index for
779 AliError("Tree with the space point arrays not initialized!");
784 Int_t nVolIds = volids->GetSize();
786 AliError("Volume IDs array is empty!");
792 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
793 UShort_t volid = (*volids)[iVolId];
795 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
797 // In case of empty index
798 if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
799 AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
802 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
806 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
811 AliTrackPointArray* array = 0;
812 fPointsTree->SetBranchAddress("SP", &array);
814 // Allocate the pointer to the space-point arrays
815 points = new AliTrackPointArray*[nArrays];
816 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
818 // Init the array used to flag already loaded tree entries
819 Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
820 for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
821 indexUsed[i] = kFALSE;
823 // Start the loop over the volume ids
825 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
826 UShort_t volid = (*volids)[iVolId];
828 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
830 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
831 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
834 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
837 Int_t entry = (*index)[iArrayId];
838 if (indexUsed[entry] == kTRUE) {
842 fPointsTree->GetEvent(entry);
844 AliWarning("Wrong space point array index!");
847 indexUsed[entry] = kTRUE;
849 // Get the space-point array
850 Int_t nPoints = array->GetNPoints();
851 points[iArray] = new AliTrackPointArray(nPoints);
852 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
853 array->GetPoint(p,iPoint);
855 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
856 // check if the layer id is valid
857 if ((layer < AliGeomManager::kFirstLayer) ||
858 (layer >= AliGeomManager::kLastLayer)) {
859 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
860 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
863 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
865 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
866 layer,modnum,AliGeomManager::LayerSize(layer)));
870 // Misalignment is introduced here
871 // Switch it off in case of real
874 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
876 misalignObj->Transform(p);
878 // End of misalignment
880 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
881 alignObj->Transform(p);
882 points[iArray]->AddPoint(iPoint,&p);
894 //______________________________________________________________________________
895 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
897 // Unload track point arrays for a given
899 for (Int_t iArray = 0; iArray < n; iArray++)
900 delete points[iArray];
904 //______________________________________________________________________________
905 AliTrackFitter *AliAlignmentTracks::CreateFitter()
907 // Check if the user has already supplied
908 // a track fitter object.
909 // If not, create a default one.
911 fTrackFitter = new AliTrackFitterRieman;
916 //______________________________________________________________________________
917 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
919 // Check if the user has already supplied
920 // a track residuals minimizer object.
921 // If not, create a default one.
923 fMinimizer = new AliTrackResidualsChi2;
928 //______________________________________________________________________________
929 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
931 // The method reads from a file a set of AliAlignObj which are
932 // then used to apply misalignments directly on the track
933 // space-points. The method is supposed to be used only for
934 // fast development and debugging of the alignment algorithms.
935 // Be careful not to use it in the case of 'real' alignment
936 // scenario since it will bias the results.
938 // Initialize the misalignment objects array
939 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
940 fMisalignObjs = new AliAlignObj**[nLayers];
941 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
942 fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
943 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
944 fMisalignObjs[iLayer][iModule] = 0x0;
947 // Open the misliagnment file and load the array with
948 // misalignment objects
949 TFile* inFile = TFile::Open(misalignObjFileName,"READ");
950 if (!inFile || !inFile->IsOpen()) {
951 AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
955 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
957 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
963 // Store the misalignment objects for further usage
964 Int_t nObjs = array->GetEntriesFast();
965 AliGeomManager::ELayerID layerId; // volume layer
966 Int_t modId; // volume ID inside the layer
967 for(Int_t i=0; i<nObjs; i++)
969 AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
970 alObj->GetVolUID(layerId,modId);
971 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
977 //________________________________________________
978 void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
981 TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
982 TClonesArray &alo=*clonesArray;
983 for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
985 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
987 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
988 new(alo[last])AliAlignObjParams(*alignObj);
992 TFile *file=new TFile(outfilename.Data(),"RECREATE");
995 alo.Write("ITSAlignObjs",TObject::kSingleKey);