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"
38 ClassImp(AliAlignmentTracks)
40 //______________________________________________________________________________
41 AliAlignmentTracks::AliAlignmentTracks():
43 fPointsFilename("AliTrackPoints.root"),
48 fIsIndexBuilt(kFALSE),
56 // Default constructor
61 //______________________________________________________________________________
62 AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
64 fPointsFilename("AliTrackPoints.root"),
69 fIsIndexBuilt(kFALSE),
77 // Constructor in the case
78 // the user provides an already
79 // built TChain with ESD trees
85 //______________________________________________________________________________
86 AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
87 fESDChain(new TChain(esdtreename)),
88 fPointsFilename("AliTrackPoints.root"),
93 fIsIndexBuilt(kFALSE),
101 // Constructor in the case
102 // the user provides a single ESD file
103 // or a directory containing ESD files
104 fESDChain->Add(esdfilename);
111 //______________________________________________________________________________
112 AliAlignmentTracks::~AliAlignmentTracks()
115 if (fESDChain) delete fESDChain;
123 if (fPointsFile) fPointsFile->Close();
126 //______________________________________________________________________________
127 void AliAlignmentTracks::AddESD(TChain *esdchain)
129 // Add a chain with ESD files
131 fESDChain->Add(esdchain);
133 fESDChain = esdchain;
136 //______________________________________________________________________________
137 void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
139 // Add a single file or
140 // a directory to the chain
141 // with the ESD files
143 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
145 fESDChain = new TChain(esdtreename);
146 fESDChain->Add(esdfilename);
151 //________________________________________________________________________
152 void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
155 Float_t minAngleWrtITSModulePlanes,
156 Float_t minMom,Float_t maxMom,
157 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
158 Float_t minSinTheta,Float_t maxSinTheta)
160 // Analyzes and filters ESD tracks
161 // Stores the selected track space points
162 // into the output file
164 if (!fESDChain) return;
166 AliESDEvent *esd = new AliESDEvent();
167 esd->ReadFromTree(fESDChain);
168 AliESDfriend *esdf = 0;
169 fESDChain->SetBranchStatus("ESDfriend*",1);
170 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
172 // Open the output file
173 if (fPointsFilename.IsNull()) {
174 AliWarning("Incorrect output filename!");
178 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
179 if (!pointsFile || !pointsFile->IsOpen()) {
180 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
184 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
185 const AliTrackPointArray *array = 0;
186 AliTrackPointArray *array2 = 0;
187 if(onlyITS) { // only ITS AliTrackPoints
188 pointsTree->Branch("SP","AliTrackPointArray", &array2);
190 pointsTree->Branch("SP","AliTrackPointArray", &array);
194 while (fESDChain->GetEntry(ievent++)) {
197 esd->SetESDfriend(esdf); //Attach the friend to the ESD
199 Int_t ntracks = esd->GetNumberOfTracks();
200 for (Int_t itrack=0; itrack < ntracks; itrack++) {
201 AliESDtrack * track = esd->GetTrack(itrack);
202 if (!track) continue;
204 if(track->GetNcls(0) < minITSpts) continue;
206 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
207 Float_t abssinphi = TMath::Abs(TMath::Sin(track->GetAlpha()+TMath::ASin(track->GetSnp())));
208 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
209 Float_t sintheta = TMath::Sin(0.5*TMath::Pi()-TMath::ATan(track->GetTgl()));
210 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
214 array = track->GetTrackPointArray();
217 Bool_t layerOK[6]={kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
218 Int_t ipt,volId,modId,layerId;
220 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
221 array->GetPoint(point,ipt);
222 volId = point.GetVolumeID();
223 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
224 if(layerId>6) continue;
225 // check minAngleWrtITSModulePlanes
227 Double_t p[3]; track->GetDirection(p);
228 TVector3 pvec(p[0],p[1],p[2]);
229 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
230 TVector3 normvec(rot[1],rot[4],rot[7]);
231 Double_t angle = pvec.Angle(normvec);
232 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
233 angle = 0.5*TMath::Pi()-angle;
234 if(angle<minAngleWrtITSModulePlanes) {
235 layerOK[layerId-1]=kFALSE;
241 if(jpt < minITSpts) continue;
242 array2 = new AliTrackPointArray(jpt);
244 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
245 array->GetPoint(point,ipt);
246 volId = point.GetVolumeID();
247 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
248 if(layerId>6 || !layerOK[layerId-1]) continue;
249 array2->AddPoint(jpt,&point);
258 if (!pointsTree->Write()) {
259 AliWarning("Can't write the tree with track point arrays!");
268 //_____________________________________________________________________________
269 void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
270 Int_t minITSpts,Float_t maxMatchingAngle,
272 Float_t minAngleWrtITSModulePlanes,
273 Float_t minMom,Float_t maxMom,
274 Float_t minAbsSinPhi,Float_t maxAbsSinPhi,
275 Float_t minSinTheta,Float_t maxSinTheta)
277 // Analyzes and filters ESD tracks
278 // Merges inward and outward tracks in one single track
279 // Stores the selected track space points
280 // into the output file
282 if (!fESDChain) return;
284 AliESDEvent *esd = new AliESDEvent();
285 esd->ReadFromTree(fESDChain);
286 AliESDfriend *esdf = 0;
287 fESDChain->SetBranchStatus("ESDfriend*",1);
288 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
290 // Open the output file
291 if (fPointsFilename.IsNull()) {
292 AliWarning("Incorrect output filename!");
296 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
297 if (!pointsFile || !pointsFile->IsOpen()) {
298 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
302 TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
303 const AliTrackPointArray *array = 0;
304 AliTrackPointArray *array2 = 0;
305 pointsTree->Branch("SP","AliTrackPointArray", &array2);
308 while (fESDChain->GetEntry(ievent++)) {
311 esd->SetESDfriend(esdf); //Attach the friend to the ESD
313 Int_t ntracks = esd->GetNumberOfTracks();
314 if(ntracks<2) continue;
315 Int_t *goodtracksArray = new Int_t[ntracks];
316 Float_t *phiArray = new Float_t[ntracks];
317 Float_t *thetaArray = new Float_t[ntracks];
319 for (Int_t itrack=0; itrack < ntracks; itrack++) {
320 AliESDtrack * track = esd->GetTrack(itrack);
321 if (!track) continue;
323 if(track->GetNcls(0) < minITSpts) continue;
324 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
325 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
327 if(track->GetP()<minMom || track->GetP()>maxMom) continue;
328 Float_t abssinphi = TMath::Abs(TMath::Sin(phi));
329 if(abssinphi<minAbsSinPhi || abssinphi>maxAbsSinPhi) continue;
330 Float_t sintheta = TMath::Sin(theta);
331 if(sintheta<minSinTheta || sintheta>maxSinTheta) continue;
333 goodtracksArray[ngt]=itrack;
335 thetaArray[ngt]=theta;
340 delete [] goodtracksArray; goodtracksArray=0;
341 delete [] phiArray; phiArray=0;
342 delete [] thetaArray; thetaArray=0;
346 // check matching of the two tracks from the muon
347 Float_t min = 10000000.;
348 Int_t good1 = -1, good2 = -1;
349 for(Int_t itr1=0; itr1<ngt-1; itr1++) {
350 for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
351 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
352 if(deltatheta>maxMatchingAngle) continue;
353 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
354 if(deltaphi>maxMatchingAngle) continue;
355 //printf("%f %f %f %f\n",deltaphi,deltatheta,thetaArray[itr1],thetaArray[itr2]);
356 if(deltatheta+deltaphi<min) {
357 min=deltatheta+deltaphi;
358 good1 = goodtracksArray[itr1];
359 good2 = goodtracksArray[itr2];
364 delete [] goodtracksArray; goodtracksArray=0;
365 delete [] phiArray; phiArray=0;
366 delete [] thetaArray; thetaArray=0;
368 if(good1<0) continue;
370 AliESDtrack * track1 = esd->GetTrack(good1);
371 AliESDtrack * track2 = esd->GetTrack(good2);
374 Int_t ipt,volId,modId,layerId;
376 Bool_t layerOK[6][2];
377 for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kTRUE;
378 array = track1->GetTrackPointArray();
379 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
380 array->GetPoint(point,ipt);
382 volId = point.GetVolumeID();
383 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
384 if(layerId>6) continue;
385 // check minAngleWrtITSModulePlanes
387 Double_t p[3]; track1->GetDirection(p);
388 TVector3 pvec(p[0],p[1],p[2]);
389 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
390 TVector3 normvec(rot[1],rot[4],rot[7]);
391 Double_t angle = pvec.Angle(normvec);
392 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
393 angle = 0.5*TMath::Pi()-angle;
394 if(angle<minAngleWrtITSModulePlanes) {
395 layerOK[layerId-1][0]=kFALSE;
402 array = track2->GetTrackPointArray();
403 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
404 array->GetPoint(point,ipt);
406 volId = point.GetVolumeID();
407 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
408 if(layerId>6) continue;
409 // check minAngleWrtITSModulePlanes
411 Double_t p[3]; track2->GetDirection(p);
412 TVector3 pvec(p[0],p[1],p[2]);
413 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
414 TVector3 normvec(rot[1],rot[4],rot[7]);
415 Double_t angle = pvec.Angle(normvec);
416 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
417 angle = 0.5*TMath::Pi()-angle;
418 if(angle<minAngleWrtITSModulePlanes) {
419 layerOK[layerId-1][0]=kFALSE;
427 if(jpt < 2*minITSpts) continue;
428 array2 = new AliTrackPointArray(jpt);
430 array = track1->GetTrackPointArray();
431 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
432 array->GetPoint(point,ipt);
434 volId = point.GetVolumeID();
435 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
436 if(layerId>6 || !layerOK[layerId-1][0]) continue;
438 array2->AddPoint(jpt,&point);
441 array = track2->GetTrackPointArray();
442 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
443 array->GetPoint(point,ipt);
445 volId = point.GetVolumeID();
446 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
447 if(layerId>6 || !layerOK[layerId-1][1]) continue;
449 array2->AddPoint(jpt,&point);
456 if (!pointsTree->Write()) {
457 AliWarning("Can't write the tree with track point arrays!");
465 //______________________________________________________________________________
466 void AliAlignmentTracks::ProcessESD(TSelector *selector)
468 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
471 //______________________________________________________________________________
472 void AliAlignmentTracks::BuildIndex()
474 // Build index of points tree entries
475 // Used for access based on the volume IDs
476 if (fIsIndexBuilt) return;
478 fIsIndexBuilt = kTRUE;
480 // Dummy object is created in order
481 // to initialize the volume paths
482 AliAlignObjParams alobj;
484 fPointsFile = TFile::Open(fPointsFilename);
485 if (!fPointsFile || !fPointsFile->IsOpen()) {
486 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
490 // AliTrackPointArray* array = new AliTrackPointArray;
491 AliTrackPointArray* array = 0;
492 fPointsTree = (TTree*) fPointsFile->Get("spTree");
494 AliWarning("No pointsTree found!");
497 fPointsTree->SetBranchAddress("SP", &array);
499 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
500 for (Int_t iArray = 0; iArray < nArrays; iArray++)
502 fPointsTree->GetEvent(iArray);
503 if (!array) continue;
504 for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
505 UShort_t volId = array->GetVolumeID()[ipoint];
506 // check if the volId is valid
507 if (!AliGeomManager::SymName(volId)) {
508 AliError(Form("The volume id %d has no default volume name !",
513 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId)
514 - AliGeomManager::kFirstLayer;
515 if (!fArrayIndex[layerId][modId]) {
516 //first entry for this volume
517 fArrayIndex[layerId][modId] = new TArrayI(1000);
520 Int_t size = fArrayIndex[layerId][modId]->GetSize();
521 // If needed allocate new size
522 if (fLastIndex[layerId][modId] >= size)
523 fArrayIndex[layerId][modId]->Set(size + 1000);
526 // Check if the index is already filled
527 Bool_t fillIndex = kTRUE;
528 if (fLastIndex[layerId][modId] != 0) {
529 if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
532 // Fill the index array and store last filled index
534 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
535 fLastIndex[layerId][modId]++;
542 //______________________________________________________________________________
543 void AliAlignmentTracks::InitIndex()
545 // Initialize the index arrays
546 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
547 fLastIndex = new Int_t*[nLayers];
548 fArrayIndex = new TArrayI**[nLayers];
549 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
550 fLastIndex[iLayer] = new Int_t[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
551 fArrayIndex[iLayer] = new TArrayI*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
552 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
553 fLastIndex[iLayer][iModule] = 0;
554 fArrayIndex[iLayer][iModule] = 0;
559 //______________________________________________________________________________
560 void AliAlignmentTracks::ResetIndex()
562 // Reset the value of the last filled index
563 // Do not realocate memory
565 fIsIndexBuilt = kFALSE;
567 for (Int_t iLayer = 0; iLayer < AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer; iLayer++) {
568 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
569 fLastIndex[iLayer][iModule] = 0;
574 //______________________________________________________________________________
575 void AliAlignmentTracks::DeleteIndex()
577 // Delete the index arrays
578 // Called by the destructor
579 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
580 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
581 if (fArrayIndex[iLayer][iModule]) {
582 delete fArrayIndex[iLayer][iModule];
583 fArrayIndex[iLayer][iModule] = 0;
586 delete [] fLastIndex[iLayer];
587 delete [] fArrayIndex[iLayer];
589 delete [] fLastIndex;
590 delete [] fArrayIndex;
593 //______________________________________________________________________________
594 Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
596 // Read alignment object from a file: update the alignobj already present with the one in the file
597 // To be replaced by a call to CDB
599 if(gSystem->AccessPathName(alignObjFileName,kFileExists)){
600 printf("Wrong AlignObjs File Name \n");
604 TFile *fRealign=TFile::Open(alignObjFileName);
605 if (!fRealign || !fRealign->IsOpen()) {
606 AliError(Form("Could not open Align Obj File file %s !",alignObjFileName));
609 printf("Getting TClonesArray \n");
610 TClonesArray *clnarray=(TClonesArray*)fRealign->Get(arrayName);
611 Int_t size=clnarray->GetSize();
614 for(Int_t ivol=0;ivol<size;ivol++){
615 AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
616 volid=a->GetVolUID();
618 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
619 if(iLayer<AliGeomManager::kFirstLayer||iLayer>AliGeomManager::kSSD2)continue;
620 printf("Updating volume: %d ,layer: %d module: %d \n",volid,iLayer,iModule);
621 *fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule] *= *a;
629 //______________________________________________________________________________
630 void AliAlignmentTracks::InitAlignObjs()
632 // Initialize the alignment objects array
633 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
634 fAlignObjs = new AliAlignObj**[nLayers];
635 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
636 fAlignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
637 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
638 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer+ AliGeomManager::kFirstLayer,iModule);
639 fAlignObjs[iLayer][iModule] = new AliAlignObjParams(AliGeomManager::SymName(volid),volid,0,0,0,0,0,0,kTRUE);
644 //______________________________________________________________________________
645 void AliAlignmentTracks::ResetAlignObjs()
647 // Reset the alignment objects array
648 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
649 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
650 fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
654 //______________________________________________________________________________
655 void AliAlignmentTracks::DeleteAlignObjs()
657 // Delete the alignment objects array
658 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
659 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
660 if (fAlignObjs[iLayer][iModule])
661 delete fAlignObjs[iLayer][iModule];
662 delete [] fAlignObjs[iLayer];
664 delete [] fAlignObjs;
668 Bool_t AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
669 AliGeomManager::ELayerID lastLayer,
670 AliGeomManager::ELayerID layerRangeMin,
671 AliGeomManager::ELayerID layerRangeMax,
674 // Align detector volumes within
675 // a given layer range
676 // (could be whole detector).
677 // Tracks are fitted only within
678 // the range defined by the user.
680 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
681 nModules += AliGeomManager::LayerSize(iLayer);
682 TArrayI volIds(nModules);
685 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++) {
686 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer); iModule++) {
687 UShort_t volId = AliGeomManager::LayerToVolUID(iLayer,iModule);
688 volIds.AddAt(volId,modnum);
693 Bool_t result = kFALSE;
694 while (iterations > 0) {
695 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
701 //______________________________________________________________________________
702 Bool_t AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
703 AliGeomManager::ELayerID layerRangeMin,
704 AliGeomManager::ELayerID layerRangeMax,
707 // Align detector volumes within
709 // Tracks are fitted only within
710 // the range defined by the user.
711 Int_t nModules = AliGeomManager::LayerSize(layer);
712 TArrayI volIds(nModules);
713 for (Int_t iModule = 0; iModule < nModules; iModule++) {
714 UShort_t volId = AliGeomManager::LayerToVolUID(layer,iModule);
715 volIds.AddAt(volId,iModule);
718 Bool_t result = kFALSE;
719 while (iterations > 0) {
720 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
726 //______________________________________________________________________________
727 Bool_t AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
730 // Align single detector volume to
732 // Tracks are fitted only within
733 // the second volume.
735 volIds.AddAt(volId,0);
736 TArrayI volIdsFit(1);
737 volIdsFit.AddAt(volIdFit,0);
739 Bool_t result = kFALSE;
740 while (iterations > 0) {
741 if (!(result = AlignVolumes(&volIds,&volIdsFit))) break;
747 //______________________________________________________________________________
748 Bool_t AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
749 AliGeomManager::ELayerID layerRangeMin,
750 AliGeomManager::ELayerID layerRangeMax,
753 // Align a set of detector volumes.
754 // Tracks are fitted only within
755 // the range defined by the user
756 // (by layerRangeMin and layerRangeMax)
757 // or within the set of volidsfit
758 // Repeat the procedure 'iterations' times
760 Int_t nVolIds = volids->GetSize();
762 AliError("Volume IDs array is empty!");
766 // Load only the tracks with at least one
767 // space point in the set of volume (volids)
769 AliTrackPointArray **points;
771 // Start the iterations
772 Bool_t result = kFALSE;
773 while (iterations > 0) {
774 Int_t nArrays = LoadPoints(volids, points,pointsdim);
775 if (nArrays == 0) return kFALSE;
777 AliTrackResiduals *minimizer = CreateMinimizer();
778 minimizer->SetNTracks(nArrays);
779 minimizer->InitAlignObj();
780 AliTrackFitter *fitter = CreateFitter();
781 for (Int_t iArray = 0; iArray < nArrays; iArray++) {
782 if (!points[iArray]) continue;
783 fitter->SetTrackPointArray(points[iArray], kFALSE);
784 if (fitter->Fit(volids,volidsfit,layerRangeMin,layerRangeMax) == kFALSE) continue;
785 AliTrackPointArray *pVolId,*pTrack;
786 fitter->GetTrackResiduals(pVolId,pTrack);
787 minimizer->AddTrackPointArrays(pVolId,pTrack);
789 if (!(result = minimizer->Minimize())) break;
791 // Update the alignment object(s)
792 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
793 UShort_t volid = (*volids)[iVolId];
795 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
796 AliAlignObj *alignObj = fAlignObjs[iLayer-AliGeomManager::kFirstLayer][iModule];
797 *alignObj *= *minimizer->GetAlignObj();
798 if(iterations==1)alignObj->Print("");
801 UnloadPoints(pointsdim, points);
808 //______________________________________________________________________________
809 Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points,Int_t &pointsdim)
811 // Load track point arrays with at least
812 // one space point in a given set of detector
813 // volumes (array volids).
814 // Use the already created tree index for
818 AliError("Tree with the space point arrays not initialized!");
823 Int_t nVolIds = volids->GetSize();
825 AliError("Volume IDs array is empty!");
831 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
832 UShort_t volid = (*volids)[iVolId];
834 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
836 // In case of empty index
837 if (fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule] == 0) {
838 AliWarning(Form("There are no space-points belonging to the volume which is to be aligned (Volume ID =%d)!",volid));
841 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
845 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
850 AliTrackPointArray* array = 0;
851 fPointsTree->SetBranchAddress("SP", &array);
853 // Allocate the pointer to the space-point arrays
855 points = new AliTrackPointArray*[nArrays];
856 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
858 // Init the array used to flag already loaded tree entries
859 Bool_t *indexUsed = new Bool_t[(UInt_t)fPointsTree->GetEntries()];
860 for (Int_t i = 0; i < fPointsTree->GetEntries(); i++)
861 indexUsed[i] = kFALSE;
863 // Start the loop over the volume ids
865 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
866 UShort_t volid = (*volids)[iVolId];
868 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
870 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
871 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
874 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
877 Int_t entry = (*index)[iArrayId];
878 if (indexUsed[entry] == kTRUE) {
882 fPointsTree->GetEvent(entry);
884 AliWarning("Wrong space point array index!");
887 indexUsed[entry] = kTRUE;
889 // Get the space-point array
890 Int_t nPoints = array->GetNPoints();
891 points[iArray] = new AliTrackPointArray(nPoints);
892 for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
893 array->GetPoint(p,iPoint);
895 AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(p.GetVolumeID(),modnum);
896 // check if the layer id is valid
897 if ((layer < AliGeomManager::kFirstLayer) ||
898 (layer >= AliGeomManager::kLastLayer)) {
899 AliError(Form("Layer index is invalid: %d (%d -> %d) !",
900 layer,AliGeomManager::kFirstLayer,AliGeomManager::kLastLayer-1));
903 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
905 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
906 layer,modnum,AliGeomManager::LayerSize(layer)));
910 // Misalignment is introduced here
911 // Switch it off in case of real
914 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
916 misalignObj->Transform(p);
918 // End of misalignment
921 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
922 UShort_t volp=p.GetVolumeID();
925 for (Int_t iVol = 0; iVol < nVolIds; iVol++) {
926 UShort_t vol = (*volids)[iVol];
928 alignObj->Transform(p,kFALSE);
934 if(!found)alignObj->Transform(p,fCovIsUsed);
935 points[iArray]->AddPoint(iPoint,&p);
947 //______________________________________________________________________________
948 void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
950 // Unload track point arrays for a given
952 for (Int_t iArray = 0; iArray < n; iArray++)
953 delete points[iArray];
957 //______________________________________________________________________________
958 AliTrackFitter *AliAlignmentTracks::CreateFitter()
960 // Check if the user has already supplied
961 // a track fitter object.
962 // If not, create a default one.
964 fTrackFitter = new AliTrackFitterRieman;
969 //______________________________________________________________________________
970 AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
972 // Check if the user has already supplied
973 // a track residuals minimizer object.
974 // If not, create a default one.
976 fMinimizer = new AliTrackResidualsChi2;
981 //______________________________________________________________________________
982 Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
984 // The method reads from a file a set of AliAlignObj which are
985 // then used to apply misalignments directly on the track
986 // space-points. The method is supposed to be used only for
987 // fast development and debugging of the alignment algorithms.
988 // Be careful not to use it in the case of 'real' alignment
989 // scenario since it will bias the results.
991 // Initialize the misalignment objects array
992 Int_t nLayers = AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer;
993 fMisalignObjs = new AliAlignObj**[nLayers];
994 for (Int_t iLayer = 0; iLayer < (AliGeomManager::kLastLayer - AliGeomManager::kFirstLayer); iLayer++) {
995 fMisalignObjs[iLayer] = new AliAlignObj*[AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer)];
996 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++)
997 fMisalignObjs[iLayer][iModule] = 0x0;
1000 // Open the misliagnment file and load the array with
1001 // misalignment objects
1002 TFile* inFile = TFile::Open(misalignObjFileName,"READ");
1003 if (!inFile || !inFile->IsOpen()) {
1004 AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
1008 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
1010 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
1016 // Store the misalignment objects for further usage
1017 Int_t nObjs = array->GetEntriesFast();
1018 AliGeomManager::ELayerID layerId; // volume layer
1019 Int_t modId; // volume ID inside the layer
1020 for(Int_t i=0; i<nObjs; i++)
1022 AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
1023 alObj->GetVolUID(layerId,modId);
1024 if(layerId<AliGeomManager::kFirstLayer) {
1025 AliWarning(Form("Alignment object is ignored: %s",alObj->GetSymName()));
1028 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
1034 //________________________________________________
1035 void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
1038 TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
1039 TClonesArray &alo=*clonesArray;
1040 for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
1042 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
1044 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
1045 new(alo[last])AliAlignObjParams(*alignObj);
1049 TFile *file=new TFile(outfilename.Data(),"RECREATE");
1052 alo.Write("ITSAlignObjs",TObject::kSingleKey);