STEER/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 #include "AliESDfriend.h"
  38
  39 ClassImp(AliAlignmentTracks)
  40
  41 //______________________________________________________________________________
  42 AliAlignmentTracks::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 //______________________________________________________________________________
  63 AliAlignmentTracks::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 //______________________________________________________________________________
  87 AliAlignmentTracks::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 //______________________________________________________________________________
 113 AliAlignmentTracks::~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 //______________________________________________________________________________
 128 void 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 //______________________________________________________________________________
 138 void 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 //________________________________________________________________________
 153 void 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 //_____________________________________________________________________________
 270 void 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 //______________________________________________________________________________
 467 void AliAlignmentTracks::ProcessESD(TSelector *selector)
 468 {
 469   AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
 470 }
 471
 472 //______________________________________________________________________________
 473 void 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 //______________________________________________________________________________
 544 void 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 //______________________________________________________________________________
 561 void 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 //______________________________________________________________________________
 576 void 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 //______________________________________________________________________________
 595 Bool_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 //______________________________________________________________________________
 631 void 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 //______________________________________________________________________________
 646 void 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 //______________________________________________________________________________
 656 void 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
 669 Bool_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 //______________________________________________________________________________
 703 Bool_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 //______________________________________________________________________________
 728 Bool_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 //______________________________________________________________________________
 749 Bool_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 //______________________________________________________________________________
 816 Int_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 //______________________________________________________________________________
 955 void 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 //______________________________________________________________________________
 965 AliTrackFitter *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 //______________________________________________________________________________
 977 AliTrackResiduals *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 //______________________________________________________________________________
 989 Bool_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 //________________________________________________
1042 void 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 }