]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - STEER/AliAlignmentTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Replaced AliInfo with AliDebug
[u/mrichter/AliRoot.git] / STEER / AliAlignmentTracks.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------
17// Implementation of the alignment steering class
18// It provides an access to the track space points
19// written along the esd tracks. The class enables
20// the user to plug any track fitter (deriving from
21// AliTrackFitter class) and minimization fo the
22// track residual sums (deriving from the AliTrackResiduals).
23//-----------------------------------------------------------------
24
25#include <TChain.h>
26#include <TFile.h>
27#include <TVector3.h>
28#include <TSystem.h>
29
30#include "AliAlignmentTracks.h"
31#include "AliTrackPointArray.h"
32#include "AliAlignObjParams.h"
33#include "AliTrackFitterRieman.h"
34#include "AliTrackResidualsChi2.h"
35#include "AliESDEvent.h"
36#include "AliLog.h"
37
38ClassImp(AliAlignmentTracks)
39
40//______________________________________________________________________________
41AliAlignmentTracks::AliAlignmentTracks():
42 fESDChain(0),
43 fPointsFilename("AliTrackPoints.root"),
44 fPointsFile(0),
45 fPointsTree(0),
46 fLastIndex(0),
47 fArrayIndex(0),
48 fIsIndexBuilt(kFALSE),
49 fAlignObjs(0),
50 fMisalignObjs(0),
51 fTrackFitter(0),
52 fMinimizer(0),
53 fDoUpdate(kTRUE),
54 fCovIsUsed(kFALSE)
55{
56 // Default constructor
57 InitIndex();
58 InitAlignObjs();
59}
60
61//______________________________________________________________________________
62AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
63 fESDChain(esdchain),
64 fPointsFilename("AliTrackPoints.root"),
65 fPointsFile(0),
66 fPointsTree(0),
67 fLastIndex(0),
68 fArrayIndex(0),
69 fIsIndexBuilt(kFALSE),
70 fAlignObjs(0),
71 fMisalignObjs(0),
72 fTrackFitter(0),
73 fMinimizer(0),
74 fDoUpdate(kTRUE),
75 fCovIsUsed(kFALSE)
76{
77 // Constructor in the case
78 // the user provides an already
79 // built TChain with ESD trees
80 InitIndex();
81 InitAlignObjs();
82}
83
84
85//______________________________________________________________________________
86AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
87 fESDChain(new TChain(esdtreename)),
88 fPointsFilename("AliTrackPoints.root"),
89 fPointsFile(0),
90 fPointsTree(0),
91 fLastIndex(0),
92 fArrayIndex(0),
93 fIsIndexBuilt(kFALSE),
94 fAlignObjs(0),
95 fMisalignObjs(0),
96 fTrackFitter(0),
97 fMinimizer(0),
98 fDoUpdate(kTRUE),
99 fCovIsUsed(kFALSE)
100{
101 // Constructor in the case
102 // the user provides a single ESD file
103 // or a directory containing ESD files
104 fESDChain->Add(esdfilename);
105
106 InitIndex();
107 InitAlignObjs();
108}
109
110
111//______________________________________________________________________________
112AliAlignmentTracks::~AliAlignmentTracks()
113{
114 // Destructor
115 if (fESDChain) delete fESDChain;
116
117 DeleteIndex();
118 DeleteAlignObjs();
119
120 delete fTrackFitter;
121 delete fMinimizer;
122
123 if (fPointsFile) fPointsFile->Close();
124}
125
126//______________________________________________________________________________
127void AliAlignmentTracks::AddESD(TChain *esdchain)
128{
129 // Add a chain with ESD files
130 if (fESDChain)
131 fESDChain->Add(esdchain);
132 else
133 fESDChain = esdchain;
134}
135
136//______________________________________________________________________________
137void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
138{
139 // Add a single file or
140 // a directory to the chain
141 // with the ESD files
142 if (fESDChain)
143 fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
144 else {
145 fESDChain = new TChain(esdtreename);
146 fESDChain->Add(esdfilename);
147 }
148}
149
150
151//________________________________________________________________________
152void AliAlignmentTracks::ProcessESD(Bool_t onlyITS,
153 Int_t minITSpts,
154 Bool_t cuts,
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)
159{
160 // Analyzes and filters ESD tracks
161 // Stores the selected track space points
162 // into the output file
163
164 if (!fESDChain) return;
165
166 AliESDEvent *esd = new AliESDEvent();
167 esd->ReadFromTree(fESDChain);
168 AliESDfriend *esdf = 0;
169 fESDChain->SetBranchStatus("ESDfriend*",1);
170 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
171
172 // Open the output file
173 if (fPointsFilename.IsNull()) {
174 AliWarning("Incorrect output filename!");
175 return;
176 }
177
178 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
179 if (!pointsFile || !pointsFile->IsOpen()) {
180 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
181 return;
182 }
183
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);
189 } else {
190 pointsTree->Branch("SP","AliTrackPointArray", &array);
191 }
192
193 Int_t ievent = 0;
194 while (fESDChain->GetEntry(ievent++)) {
195 if (!esd) break;
196
197 esd->SetESDfriend(esdf); //Attach the friend to the ESD
198
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;
203
204 if(track->GetNcls(0) < minITSpts) continue;
205 if(cuts) {
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;
211 }
212
213 AliTrackPoint point;
214 array = track->GetTrackPointArray();
215
216 if(onlyITS) {
217 Bool_t layerOK[6]={kTRUE,kTRUE,kTRUE,kTRUE,kTRUE,kTRUE};
218 Int_t ipt,volId,modId,layerId;
219 Int_t jpt=0;
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
226 if(cuts) {
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;
236 continue;
237 }
238 }
239 jpt++;
240 }
241 if(jpt < minITSpts) continue;
242 array2 = new AliTrackPointArray(jpt);
243 jpt=0;
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);
250 jpt++;
251 }
252 } // end if(onlyITS)
253
254 pointsTree->Fill();
255 }
256 }
257
258 if (!pointsTree->Write()) {
259 AliWarning("Can't write the tree with track point arrays!");
260 return;
261 }
262
263 pointsFile->Close();
264
265 return;
266}
267
268//_____________________________________________________________________________
269void AliAlignmentTracks::ProcessESDCosmics(Bool_t onlyITS,
270 Int_t minITSpts,Float_t maxMatchingAngle,
271 Bool_t cuts,
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)
276{
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
281
282 if (!fESDChain) return;
283
284 AliESDEvent *esd = new AliESDEvent();
285 esd->ReadFromTree(fESDChain);
286 AliESDfriend *esdf = 0;
287 fESDChain->SetBranchStatus("ESDfriend*",1);
288 fESDChain->SetBranchAddress("ESDfriend.",&esdf);
289
290 // Open the output file
291 if (fPointsFilename.IsNull()) {
292 AliWarning("Incorrect output filename!");
293 return;
294 }
295
296 TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
297 if (!pointsFile || !pointsFile->IsOpen()) {
298 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
299 return;
300 }
301
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);
306
307 Int_t ievent = 0;
308 while (fESDChain->GetEntry(ievent++)) {
309 if (!esd) break;
310
311 esd->SetESDfriend(esdf); //Attach the friend to the ESD
312
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];
318 Int_t ngt=0;
319 for (Int_t itrack=0; itrack < ntracks; itrack++) {
320 AliESDtrack * track = esd->GetTrack(itrack);
321 if (!track) continue;
322
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());
326 if(cuts) {
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;
332 }
333 goodtracksArray[ngt]=itrack;
334 phiArray[ngt]=phi;
335 thetaArray[ngt]=theta;
336 ngt++;
337 }
338
339 if(ngt<2) {
340 delete [] goodtracksArray; goodtracksArray=0;
341 delete [] phiArray; phiArray=0;
342 delete [] thetaArray; thetaArray=0;
343 continue;
344 }
345
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];
360 }
361 }
362 }
363
364 delete [] goodtracksArray; goodtracksArray=0;
365 delete [] phiArray; phiArray=0;
366 delete [] thetaArray; thetaArray=0;
367
368 if(good1<0) continue;
369
370 AliESDtrack * track1 = esd->GetTrack(good1);
371 AliESDtrack * track2 = esd->GetTrack(good2);
372
373 AliTrackPoint point;
374 Int_t ipt,volId,modId,layerId;
375 Int_t jpt=0;
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);
381 if(onlyITS) {
382 volId = point.GetVolumeID();
383 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
384 if(layerId>6) continue;
385 // check minAngleWrtITSModulePlanes
386 if(cuts) {
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;
396 continue;
397 }
398 }
399 }
400 jpt++;
401 }
402 array = track2->GetTrackPointArray();
403 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
404 array->GetPoint(point,ipt);
405 if(onlyITS) {
406 volId = point.GetVolumeID();
407 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
408 if(layerId>6) continue;
409 // check minAngleWrtITSModulePlanes
410 if(cuts) {
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;
420 continue;
421 }
422 }
423 }
424 jpt++;
425 }
426
427 if(jpt < 2*minITSpts) continue;
428 array2 = new AliTrackPointArray(jpt);
429 jpt=0;
430 array = track1->GetTrackPointArray();
431 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
432 array->GetPoint(point,ipt);
433 if(onlyITS) {
434 volId = point.GetVolumeID();
435 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
436 if(layerId>6 || !layerOK[layerId-1][0]) continue;
437 }
438 array2->AddPoint(jpt,&point);
439 jpt++;
440 }
441 array = track2->GetTrackPointArray();
442 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
443 array->GetPoint(point,ipt);
444 if(onlyITS) {
445 volId = point.GetVolumeID();
446 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
447 if(layerId>6 || !layerOK[layerId-1][1]) continue;
448 }
449 array2->AddPoint(jpt,&point);
450 jpt++;
451 }
452
453 pointsTree->Fill();
454 }
455
456 if (!pointsTree->Write()) {
457 AliWarning("Can't write the tree with track point arrays!");
458 return;
459 }
460
461 pointsFile->Close();
462 return;
463}
464
465//______________________________________________________________________________
466void AliAlignmentTracks::ProcessESD(TSelector *selector)
467{
468 AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
469}
470
471//______________________________________________________________________________
472void AliAlignmentTracks::BuildIndex()
473{
474 // Build index of points tree entries
475 // Used for access based on the volume IDs
476 if (fIsIndexBuilt) return;
477
478 fIsIndexBuilt = kTRUE;
479
480 // Dummy object is created in order
481 // to initialize the volume paths
482 AliAlignObjParams alobj;
483
484 fPointsFile = TFile::Open(fPointsFilename);
485 if (!fPointsFile || !fPointsFile->IsOpen()) {
486 AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
487 return;
488 }
489
490 // AliTrackPointArray* array = new AliTrackPointArray;
491 AliTrackPointArray* array = 0;
492 fPointsTree = (TTree*) fPointsFile->Get("spTree");
493 if (!fPointsTree) {
494 AliWarning("No pointsTree found!");
495 return;
496 }
497 fPointsTree->SetBranchAddress("SP", &array);
498
499 Int_t nArrays = (Int_t)fPointsTree->GetEntries();
500 for (Int_t iArray = 0; iArray < nArrays; iArray++)
501 {
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 !",
509 volId));
510 continue;
511 }
512 Int_t modId;
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);
518 }
519 else {
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);
524 }
525
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)
530 fillIndex = kFALSE;
531 }
532 // Fill the index array and store last filled index
533 if (fillIndex) {
534 (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
535 fLastIndex[layerId][modId]++;
536 }
537 }
538 }
539
540}
541
542//______________________________________________________________________________
543void AliAlignmentTracks::InitIndex()
544{
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;
555 }
556 }
557}
558
559//______________________________________________________________________________
560void AliAlignmentTracks::ResetIndex()
561{
562 // Reset the value of the last filled index
563 // Do not realocate memory
564
565 fIsIndexBuilt = kFALSE;
566
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;
570 }
571 }
572}
573
574//______________________________________________________________________________
575void AliAlignmentTracks::DeleteIndex()
576{
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;
584 }
585 }
586 delete [] fLastIndex[iLayer];
587 delete [] fArrayIndex[iLayer];
588 }
589 delete [] fLastIndex;
590 delete [] fArrayIndex;
591}
592
593//______________________________________________________________________________
594Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
595{
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
598
599 if(gSystem->AccessPathName(alignObjFileName,kFileExists)){
600 printf("Wrong AlignObjs File Name \n");
601 return kFALSE;
602 }
603
604 TFile *fRealign=TFile::Open(alignObjFileName);
605 if (!fRealign || !fRealign->IsOpen()) {
606 AliError(Form("Could not open Align Obj File file %s !",alignObjFileName));
607 return kFALSE;
608 }
609 printf("Getting TClonesArray \n");
610 TClonesArray *clnarray=(TClonesArray*)fRealign->Get(arrayName);
611 Int_t size=clnarray->GetSize();
612 UShort_t volid;
613
614 for(Int_t ivol=0;ivol<size;ivol++){
615 AliAlignObjParams *a=(AliAlignObjParams*)clnarray->At(ivol);
616 volid=a->GetVolUID();
617 Int_t iModule;
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;
622 }
623
624 delete clnarray;
625 fRealign->Close();
626 return kTRUE;
627}
628
629//______________________________________________________________________________
630void AliAlignmentTracks::InitAlignObjs()
631{
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);
640 }
641 }
642}
643
644//______________________________________________________________________________
645void AliAlignmentTracks::ResetAlignObjs()
646{
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);
651 }
652}
653
654//______________________________________________________________________________
655void AliAlignmentTracks::DeleteAlignObjs()
656{
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];
663 }
664 delete [] fAlignObjs;
665 fAlignObjs = 0;
666}
667
668Bool_t AliAlignmentTracks::AlignDetector(AliGeomManager::ELayerID firstLayer,
669 AliGeomManager::ELayerID lastLayer,
670 AliGeomManager::ELayerID layerRangeMin,
671 AliGeomManager::ELayerID layerRangeMax,
672 Int_t iterations)
673{
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.
679 Int_t nModules = 0;
680 for (Int_t iLayer = firstLayer; iLayer <= lastLayer; iLayer++)
681 nModules += AliGeomManager::LayerSize(iLayer);
682 TArrayI volIds(nModules);
683
684 Int_t modnum = 0;
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);
689 modnum++;
690 }
691 }
692
693 Bool_t result = kFALSE;
694 while (iterations > 0) {
695 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
696 iterations--;
697 }
698 return result;
699}
700
701//______________________________________________________________________________
702Bool_t AliAlignmentTracks::AlignLayer(AliGeomManager::ELayerID layer,
703 AliGeomManager::ELayerID layerRangeMin,
704 AliGeomManager::ELayerID layerRangeMax,
705 Int_t iterations)
706{
707 // Align detector volumes within
708 // a given layer.
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);
716 }
717
718 Bool_t result = kFALSE;
719 while (iterations > 0) {
720 if (!(result = AlignVolumes(&volIds,0x0,layerRangeMin,layerRangeMax))) break;
721 iterations--;
722 }
723 return result;
724}
725
726//______________________________________________________________________________
727Bool_t AliAlignmentTracks::AlignVolume(UShort_t volId, UShort_t volIdFit,
728 Int_t iterations)
729{
730 // Align single detector volume to
731 // another volume.
732 // Tracks are fitted only within
733 // the second volume.
734 TArrayI volIds(1);
735 volIds.AddAt(volId,0);
736 TArrayI volIdsFit(1);
737 volIdsFit.AddAt(volIdFit,0);
738
739 Bool_t result = kFALSE;
740 while (iterations > 0) {
741 if (!(result = AlignVolumes(&volIds,&volIdsFit))) break;
742 iterations--;
743 }
744 return result;
745}
746
747//______________________________________________________________________________
748Bool_t AliAlignmentTracks::AlignVolumes(const TArrayI *volids, const TArrayI *volidsfit,
749 AliGeomManager::ELayerID layerRangeMin,
750 AliGeomManager::ELayerID layerRangeMax,
751 Int_t iterations)
752{
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
759
760 Int_t nVolIds = volids->GetSize();
761 if (nVolIds == 0) {
762 AliError("Volume IDs array is empty!");
763 return kFALSE;
764 }
765
766 // Load only the tracks with at least one
767 // space point in the set of volume (volids)
768 BuildIndex();
769 AliTrackPointArray **points;
770 Int_t pointsdim;
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;
776
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);
788 }
789 if (!(result = minimizer->Minimize())) break;
790
791 // Update the alignment object(s)
792 if (fDoUpdate) for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
793 UShort_t volid = (*volids)[iVolId];
794 Int_t iModule;
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("");
799 }
800
801 UnloadPoints(pointsdim, points);
802
803 iterations--;
804 }
805 return result;
806}
807
808//______________________________________________________________________________
809Int_t AliAlignmentTracks::LoadPoints(const TArrayI *volids, AliTrackPointArray** &points,Int_t &pointsdim)
810{
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
815 // fast access.
816
817 if (!fPointsTree) {
818 AliError("Tree with the space point arrays not initialized!");
819 points = 0;
820 return 0;
821 }
822
823 Int_t nVolIds = volids->GetSize();
824 if (nVolIds == 0) {
825 AliError("Volume IDs array is empty!");
826 points = 0;
827 return 0;
828 }
829
830 Int_t nArrays = 0;
831 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
832 UShort_t volid = (*volids)[iVolId];
833 Int_t iModule;
834 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
835
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));
839 continue;
840 }
841 nArrays += fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
842 }
843
844 if (nArrays == 0) {
845 AliError("There are no space-points belonging to all of the volumes which are to be aligned!");
846 points = 0x0;
847 return 0;
848 }
849
850 AliTrackPointArray* array = 0;
851 fPointsTree->SetBranchAddress("SP", &array);
852
853 // Allocate the pointer to the space-point arrays
854 pointsdim=nArrays;
855 points = new AliTrackPointArray*[nArrays];
856 for (Int_t i = 0; i < nArrays; i++) points[i] = 0x0;
857
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;
862
863 // Start the loop over the volume ids
864 Int_t iArray = 0;
865 for (Int_t iVolId = 0; iVolId < nVolIds; iVolId++) {
866 UShort_t volid = (*volids)[iVolId];
867 Int_t iModule;
868 AliGeomManager::ELayerID iLayer = AliGeomManager::VolUIDToLayer(volid,iModule);
869
870 Int_t nArraysId = fLastIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
871 TArrayI *index = fArrayIndex[iLayer-AliGeomManager::kFirstLayer][iModule];
872 AliTrackPoint p;
873
874 for (Int_t iArrayId = 0; iArrayId < nArraysId; iArrayId++) {
875
876 // Get tree entry
877 Int_t entry = (*index)[iArrayId];
878 if (indexUsed[entry] == kTRUE) {
879 nArrays--;
880 continue;
881 }
882 fPointsTree->GetEvent(entry);
883 if (!array) {
884 AliWarning("Wrong space point array index!");
885 continue;
886 }
887 indexUsed[entry] = kTRUE;
888
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);
894 Int_t modnum;
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));
901 continue;
902 }
903 if ((modnum >= AliGeomManager::LayerSize(layer)) ||
904 (modnum < 0)) {
905 AliError(Form("Module number inside layer %d is invalid: %d (0 -> %d)",
906 layer,modnum,AliGeomManager::LayerSize(layer)));
907 continue;
908 }
909
910 // Misalignment is introduced here
911 // Switch it off in case of real
912 // alignment job!
913 if (fMisalignObjs) {
914 AliAlignObj *misalignObj = fMisalignObjs[layer-AliGeomManager::kFirstLayer][modnum];
915 if (misalignObj)
916 misalignObj->Transform(p);
917 }
918 // End of misalignment
919
920
921 AliAlignObj *alignObj = fAlignObjs[layer-AliGeomManager::kFirstLayer][modnum];
922 UShort_t volp=p.GetVolumeID();
923 Bool_t found=kFALSE;
924 if(fCovIsUsed){
925 for (Int_t iVol = 0; iVol < nVolIds; iVol++) {
926 UShort_t vol = (*volids)[iVol];
927 if(volp==vol){
928 alignObj->Transform(p,kFALSE);
929 found=kTRUE;
930 break;
931 }
932 }
933 }
934 if(!found)alignObj->Transform(p,fCovIsUsed);
935 points[iArray]->AddPoint(iPoint,&p);
936 }
937 iArray++;
938 }
939 }
940
941
942 delete [] indexUsed;
943
944 return nArrays;
945}
946
947//______________________________________________________________________________
948void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
949{
950 // Unload track point arrays for a given
951 // detector volume
952 for (Int_t iArray = 0; iArray < n; iArray++)
953 delete points[iArray];
954 delete [] points;
955}
956
957//______________________________________________________________________________
958AliTrackFitter *AliAlignmentTracks::CreateFitter()
959{
960 // Check if the user has already supplied
961 // a track fitter object.
962 // If not, create a default one.
963 if (!fTrackFitter)
964 fTrackFitter = new AliTrackFitterRieman;
965
966 return fTrackFitter;
967}
968
969//______________________________________________________________________________
970AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
971{
972 // Check if the user has already supplied
973 // a track residuals minimizer object.
974 // If not, create a default one.
975 if (!fMinimizer)
976 fMinimizer = new AliTrackResidualsChi2;
977
978 return fMinimizer;
979}
980
981//______________________________________________________________________________
982Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
983{
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.
990
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;
998 }
999
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));
1005 return kFALSE;
1006 }
1007
1008 TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
1009 if (!array) {
1010 AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
1011 inFile->Close();
1012 return kFALSE;
1013 }
1014 inFile->Close();
1015
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++)
1021 {
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()));
1026 continue;
1027 }
1028 fMisalignObjs[layerId-AliGeomManager::kFirstLayer][modId] = alObj;
1029 }
1030 return kTRUE;
1031}
1032
1033
1034//________________________________________________
1035void AliAlignmentTracks::WriteRealignObjArray(TString outfilename,AliGeomManager::ELayerID layerRangeMin,AliGeomManager::ELayerID layerRangeMax){
1036
1037 Int_t last=0;
1038 TClonesArray *clonesArray=new TClonesArray("AliAlignObjParams",2200);
1039 TClonesArray &alo=*clonesArray;
1040 for (Int_t iLayer = layerRangeMin-AliGeomManager::kFirstLayer; iLayer <= (layerRangeMax - AliGeomManager::kFirstLayer);iLayer++) {
1041
1042 for (Int_t iModule = 0; iModule < AliGeomManager::LayerSize(iLayer + AliGeomManager::kFirstLayer); iModule++) {
1043
1044 AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
1045 new(alo[last])AliAlignObjParams(*alignObj);
1046 last++;
1047 }
1048 }
1049 TFile *file=new TFile(outfilename.Data(),"RECREATE");
1050 file->cd();
1051
1052 alo.Write("ITSAlignObjs",TObject::kSingleKey);
1053 file->Close();
1054
1055
1056 delete clonesArray;
1057 return;
1058}