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 //-----------------------------------------------------------------
30 #include "AliAlignmentTracks.h"
31 #include "AliTrackPointArray.h"
32 #include "AliAlignObjParams.h"
33 #include "AliTrackFitterRieman.h"
34 #include "AliTrackResidualsChi2.h"
35 #include "AliESDEvent.h"
37 #include "AliESDfriend.h"
39 ClassImp(AliAlignmentTracks)
41 //______________________________________________________________________________
42 AliAlignmentTracks::AliAlignmentTracks():
44 fPointsFilename("AliTrackPoints.root"),
49 fIsIndexBuilt(kFALSE),
57 // Default constructor
62 //______________________________________________________________________________
63 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
65 fPointsFilename("AliTrackPoints.root"),
70 fIsIndexBuilt(kFALSE),
78 // Constructor in the case
79 // the user provides an already
80 // built TChain with ESD trees
86 //______________________________________________________________________________
87 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
88 fESDChain(new TChain(esdtreename)),
89 fPointsFilename("AliTrackPoints.root"),
94 fIsIndexBuilt(kFALSE),
102 // Constructor in the case
103 // the user provides a single ESD file
104 // or a directory containing ESD files
105 fESDChain->Add(esdfilename);
112 //______________________________________________________________________________
113 AliAlignmentTracks::~AliAlignmentTracks()
116 if (fESDChain) delete fESDChain;
124 if (fPointsFile) fPointsFile->Close();
127 //______________________________________________________________________________
128 void AliAlignmentTracks::AddESD(TChain *esdchain)
130 // Add a chain with ESD files
132 fESDChain->Add(esdchain);
134 fESDChain = esdchain;
137 //______________________________________________________________________________
138 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
140 // Add a single file or
141 // a directory to the chain
142 // with the ESD files
144 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
146 fESDChain = new TChain(esdtreename);
147 fESDChain->Add(esdfilename);
152 //________________________________________________________________________
153 void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
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)
161 // Analyzes and filters ESD tracks
162 // Stores the selected track space points
163 // into the output file
165 if (!fESDChain) return;
167 AliESDEvent *esd = new AliESDEvent();
168 esd->ReadFromTree(fESDChain);
169 AliESDfriend *esdf = 0;
170 fESDChain->SetBranchStatus("ESDfriend*",1);
171 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
173 // Open the output file
174 if (fPointsFilename.IsNull()) {
175 AliWarning("Incorrect output filename!");
179 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
180 if (!pointsFile || !pointsFile->IsOpen()) {
181 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
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);
191 pointsTree->Branch("SP","AliTrackPointArray", &array);
195 while (fESDChain->GetEntry(ievent++)) {
198 esd->SetESDfriend(esdf); //Attach the friend to the ESD
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;
205 if(track->GetNcls(0) < minITSpts) continue;
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;
215 array = track->GetTrackPointArray();
218 Bool_t layerOK[6]={kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
219 Int_t ipt,volId,modId,layerId;
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
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;
242 if(jpt < minITSpts) continue;
243 array2 = new AliTrackPointArray(jpt);
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);
259 if (!pointsTree->Write()) {
260 AliWarning("Can't write the tree with track point arrays!");
269 //_____________________________________________________________________________
270 void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
271 Int_t minITSpts,Float_t maxMatchingAngle,
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)
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
283 if (!fESDChain) return;
285 AliESDEvent *esd = new AliESDEvent();
286 esd->ReadFromTree(fESDChain);
287 AliESDfriend *esdf = 0;
288 fESDChain->SetBranchStatus("ESDfriend*",1);
289 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
291 // Open the output file
292 if (fPointsFilename.IsNull()) {
293 AliWarning("Incorrect output filename!");
297 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
298 if (!pointsFile || !pointsFile->IsOpen()) {
299 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
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);
309 while (fESDChain->GetEntry(ievent++)) {
312 esd->SetESDfriend(esdf); //Attach the friend to the ESD
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];
320 for (Int_t itrack=0; itrack < ntracks; itrack++) {
321 AliESDtrack * track = esd->GetTrack(itrack);
322 if (!track) continue;
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());
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;
334 goodtracksArray[ngt]=itrack;
336 thetaArray[ngt]=theta;
341 delete [] goodtracksArray; goodtracksArray=0;
342 delete [] phiArray; phiArray=0;
343 delete [] thetaArray; thetaArray=0;
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];
365 delete [] goodtracksArray; goodtracksArray=0;
366 delete [] phiArray; phiArray=0;
367 delete [] thetaArray; thetaArray=0;
369 if(good1<0) continue;
371 AliESDtrack * track1 = esd->GetTrack(good1);
372 AliESDtrack * track2 = esd->GetTrack(good2);
375 Int_t ipt,volId,modId,layerId;
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);
383 volId = point.GetVolumeID();
384 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
385 if(layerId>6) continue;
386 // check minAngleWrtITSModulePlanes
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;
403 array = track2->GetTrackPointArray();
404 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
405 array->GetPoint(point,ipt);
407 volId = point.GetVolumeID();
408 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
409 if(layerId>6) continue;
410 // check minAngleWrtITSModulePlanes
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;
428 if(jpt < 2*minITSpts) continue;
429 array2 = new AliTrackPointArray(jpt);
431 array = track1->GetTrackPointArray();
432 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
433 array->GetPoint(point,ipt);
435 volId = point.GetVolumeID();
436 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
437 if(layerId>6 || !layerOK[layerId-1][0]) continue;
439 array2->AddPoint(jpt,&point);
442 array = track2->GetTrackPointArray();
443 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
444 array->GetPoint(point,ipt);
446 volId = point.GetVolumeID();
447 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
448 if(layerId>6 || !layerOK[layerId-1][1]) continue;
450 array2->AddPoint(jpt,&point);
457 if (!pointsTree->Write()) {
458 AliWarning("Can't write the tree with track point arrays!");
466 //______________________________________________________________________________
467 void AliAlignmentTracks::ProcessESD(TSelector *selector)
469 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
472 //______________________________________________________________________________
473 void AliAlignmentTracks::BuildIndex()
475 // Build index of points tree entries
476 // Used for access based on the volume IDs
477 if (fIsIndexBuilt) return;
479 fIsIndexBuilt = kTRUE;
481 // Dummy object is created in order
482 // to initialize the volume paths
483 AliAlignObjParams alobj;
485 fPointsFile = TFile::Open(fPointsFilename);
486 if (!fPointsFile || !fPointsFile->IsOpen()) {
487 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
491 // AliTrackPointArray* array = new AliTrackPointArray;
492 AliTrackPointArray* array = 0;
493 fPointsTree = (TTree*) fPointsFile->Get("spTree");
495 AliWarning("No pointsTree found!");
498 fPointsTree->SetBranchAddress("SP", &array);
500 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
501 for (Int_t iArray = 0; iArray < nArrays; iArray++)
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 !",
514 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
515 - AliGeomManager::kFirstLayer;
516 if (!fArrayIndex[layerId][modId]) {
517 //first entry for this volume
518 fArrayIndex[layerId][modId] = new TArrayI(1000);
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);
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)
533 // Fill the index array and store last filled index
535 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
536 fLastIndex[layerId][modId]++;
543 //______________________________________________________________________________
544 void AliAlignmentTracks::InitIndex()
546 // Initialize the index arrays
547 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
548 fLastIndex = new Int_t*[nLayers];
549 fArrayIndex = new TArrayI**[nLayers];
550 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
551 fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
552 fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
553 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
554 fLastIndex[iLayer][iModule] = 0;
555 fArrayIndex[iLayer][iModule] = 0;
560 //______________________________________________________________________________
561 void AliAlignmentTracks::ResetIndex()
563 // Reset the value of the last filled index
564 // Do not realocate memory
566 fIsIndexBuilt = kFALSE;
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;
575 //______________________________________________________________________________
576 void AliAlignmentTracks::DeleteIndex()
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;
587 delete [] fLastIndex[iLayer];
588 delete [] fArrayIndex[iLayer];
590 delete [] fLastIndex;
591 delete [] fArrayIndex;
594 //______________________________________________________________________________
595 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
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
600 if(gSystem->AccessPathName(alignObjFileName,kFileExists)){
601 printf("Wrong AlignObjs File Name \n");
605 TFile *fRealign=TFile::Open(alignObjFileName);
606 if (!fRealign || !fRealign->IsOpen()) {
607 AliError(Form("Could not open Align Obj File file %s !",alignObjFileName));
610 printf("Getting TClonesArray \n");
611 TClonesArray *clnarray=(TClonesArray*)fRealign->Get(arrayName);
612 Int_t size=clnarray->GetSize();
615 for(Int_t ivol=0;ivol<size;ivol++){
616 AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
617 volid=a->GetVolUID();
619 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
620 if(iLayer<AliGeomManager::kFirstLayer||iLayer>AliGeomManager::kSSD2)continue;
621 printf("Updating volume: %d ,layer: %d module: %d \n",volid,iLayer,iModule);
622 *fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule] *= *a;
630 //______________________________________________________________________________
631 void AliAlignmentTracks::InitAlignObjs()
633 // Initialize the alignment objects array
634 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
635 fAlignObjs = new AliAlignObj**[nLayers];
636 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
637 fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
638 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
639 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
640 fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
645 //______________________________________________________________________________
646 void AliAlignmentTracks::ResetAlignObjs()
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);
655 //______________________________________________________________________________
656 void AliAlignmentTracks::DeleteAlignObjs()
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];
665 delete [] fAlignObjs;
669 Bool_t AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
670 AliGeomManager::ELayerID lastLayer,
671 AliGeomManager::ELayerID layerRangeMin,
672 AliGeomManager::ELayerID layerRangeMax,
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.
681 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
682 nModules += AliGeomManager::LayerSize(iLayer);
683 TArrayI volIds(nModules);
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);
694 Bool_t result = kFALSE;
695 while (iterations > 0) {
696 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
702 //______________________________________________________________________________
703 Bool_t AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
704 AliGeomManager::ELayerID layerRangeMin,
705 AliGeomManager::ELayerID layerRangeMax,
708 // Align detector volumes within
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);
719 Bool_t result = kFALSE;
720 while (iterations > 0) {
721 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
727 //______________________________________________________________________________
728 Bool_t AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
731 // Align single detector volume to
733 // Tracks are fitted only within
734 // the second volume.
736 volIds.AddAt(volId,0);
737 TArrayI volIdsFit(1);
738 volIdsFit.AddAt(volIdFit,0);
740 Bool_t result = kFALSE;
741 while (iterations > 0) {
742 if (!(result = AlignVolumes(&volIds,&volIdsFit))) break;
748 //______________________________________________________________________________
749 Bool_t AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
750 AliGeomManager::ELayerID layerRangeMin,
751 AliGeomManager::ELayerID layerRangeMax,
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
761 Int_t nVolIds = volids->GetSize();
763 AliError("Volume IDs array is empty!");
767 // Load only the tracks with at least one
768 // space point in the set of volume (volids)
770 AliTrackPointArray **points;
772 // Start the iterations
773 Bool_t result = kFALSE;
774 while (iterations > 0) {
775 Int_t nArrays = LoadPoints(volids, points,pointsdim);
777 UnloadPoints(pointsdim, points);
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);
793 if (!(result = minimizer->Minimize())) {
794 UnloadPoints(pointsdim, points);
798 // Update the alignment object(s)
799 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
800 UShort_t volid = (*volids)[iVolId];
802 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
803 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
804 *alignObj *= *minimizer->GetAlignObj();
805 if(iterations==1)alignObj->Print("");
808 UnloadPoints(pointsdim, points);
815 //______________________________________________________________________________
816 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points,Int_t &pointsdim)
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
825 AliError("Tree with the space point arrays not initialized!");
830 Int_t nVolIds = volids->GetSize();
832 AliError("Volume IDs array is empty!");
838 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
839 UShort_t volid = (*volids)[iVolId];
841 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
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));
848 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
852 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
857 AliTrackPointArray* array = 0;
858 fPointsTree->SetBranchAddress("SP", &array);
860 // Allocate the pointer to the space-point arrays
862 points = new AliTrackPointArray*[nArrays];
863 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
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;
870 // Start the loop over the volume ids
872 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
873 UShort_t volid = (*volids)[iVolId];
875 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
877 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
878 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
881 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
884 Int_t entry = (*index)[iArrayId];
885 if (indexUsed[entry] == kTRUE) {
889 fPointsTree->GetEvent(entry);
891 AliWarning("Wrong space point array index!");
894 indexUsed[entry] = kTRUE;
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);
902 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
903 // check if the layer id is valid
904 if ((layer < AliGeomManager::kFirstLayer) ||
905 (layer >= AliGeomManager::kLastLayer)) {
906 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
907 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
910 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
912 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
913 layer,modnum,AliGeomManager::LayerSize(layer)));
917 // Misalignment is introduced here
918 // Switch it off in case of real
921 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
923 misalignObj->Transform(p);
925 // End of misalignment
928 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
929 UShort_t volp=p.GetVolumeID();
932 for (Int_t iVol = 0; iVol < nVolIds; iVol++) {
933 UShort_t vol = (*volids)[iVol];
935 alignObj->Transform(p,kFALSE);
941 if(!found)alignObj->Transform(p,fCovIsUsed);
942 points[iArray]->AddPoint(iPoint,&p);
954 //______________________________________________________________________________
955 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
957 // Unload track point arrays for a given
959 for (Int_t iArray = 0; iArray < n; iArray++)
960 delete points[iArray];
964 //______________________________________________________________________________
965 AliTrackFitter *AliAlignmentTracks::CreateFitter()
967 // Check if the user has already supplied
968 // a track fitter object.
969 // If not, create a default one.
971 fTrackFitter = new AliTrackFitterRieman;
976 //______________________________________________________________________________
977 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
979 // Check if the user has already supplied
980 // a track residuals minimizer object.
981 // If not, create a default one.
983 fMinimizer = new AliTrackResidualsChi2;
988 //______________________________________________________________________________
989 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
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.
998 // Initialize the misalignment objects array
999 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
1000 fMisalignObjs = new AliAlignObj**[nLayers];
1001 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
1002 fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
1003 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
1004 fMisalignObjs[iLayer][iModule] = 0x0;
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));
1015 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
1017 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
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++)
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()));
1035 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
1041 //________________________________________________
1042 void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
1045 TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
1046 TClonesArray &alo=*clonesArray;
1047 for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
1049 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
1051 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
1052 new(alo[last])AliAlignObjParams(*alignObj);
1056 TFile *file=new TFile(outfilename.Data(),"RECREATE");
1059 alo.Write("ITSAlignObjs",TObject::kSingleKey);