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