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