STEER/AliAlignmentTracks.cxx

   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
  38 ClassImp(AliAlignmentTracks)
  39
  40 //______________________________________________________________________________
  41 AliAlignmentTracks::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 //______________________________________________________________________________
  62 AliAlignmentTracks::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 //______________________________________________________________________________
  86 AliAlignmentTracks::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 //______________________________________________________________________________
 112 AliAlignmentTracks::~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 //______________________________________________________________________________
 127 void 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 //______________________________________________________________________________
 137 void 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 //________________________________________________________________________
 152 void 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 //_____________________________________________________________________________
 269 void 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 //______________________________________________________________________________
 466 void AliAlignmentTracks::ProcessESD(TSelector *selector)
 467 {
 468   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
 469 }
 470
 471 //______________________________________________________________________________
 472 void 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 //______________________________________________________________________________
 543 void 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 //______________________________________________________________________________
 560 void 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 //______________________________________________________________________________
 575 void 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 //______________________________________________________________________________
 594 Bool_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 //______________________________________________________________________________
 630 void 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 //______________________________________________________________________________
 645 void 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 //______________________________________________________________________________
 655 void 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
 668 Bool_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 //______________________________________________________________________________
 702 Bool_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 //______________________________________________________________________________
 727 Bool_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 //______________________________________________________________________________
 748 Bool_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 //______________________________________________________________________________
 809 Int_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 //______________________________________________________________________________
 948 void 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 //______________________________________________________________________________
 958 AliTrackFitter *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 //______________________________________________________________________________
 970 AliTrackResiduals *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 //______________________________________________________________________________
 982 Bool_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 //________________________________________________
1035 void 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 }