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