[u/mrichter/AliRoot.git] / STEER / AliAlignmentTracks.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------
//   Implementation of the alignment steering class
//   It provides an access to the track space points
//   written along the esd tracks. The class enables
//   the user to plug any track fitter (deriving from
//   AliTrackFitter class) and minimization fo the
//   track residual sums (deriving from the AliTrackResiduals).
//-----------------------------------------------------------------

#include <TChain.h>

#include "AliAlignmentTracks.h"
#include "AliTrackPointArray.h"
#include "AliAlignObjAngles.h"
#include "AliTrackFitterRieman.h"
#include "AliTrackResidualsChi2.h"
#include "AliESD.h"
#include "AliLog.h"

ClassImp(AliAlignmentTracks)

//______________________________________________________________________________
AliAlignmentTracks::AliAlignmentTracks():
  fESDChain(0),
  fPointsFilename("AliTrackPoints.root"),
  fPointsFile(0),
  fPointsTree(0),
  fLastIndex(0),
  fArrayIndex(0),
  fIsIndexBuilt(kFALSE),
  fTrackFitter(0),
  fMinimizer(0)
{
  // Default constructor
  InitIndex();
  InitAlignObjs();
}

//______________________________________________________________________________
AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
  fESDChain(esdchain),
  fPointsFilename("AliTrackPoints.root"),
  fPointsFile(0),
  fPointsTree(0),
  fLastIndex(0),
  fArrayIndex(0),
  fIsIndexBuilt(kFALSE),
  fTrackFitter(0),
  fMinimizer(0)
{
  // Constructor in the case
  // the user provides an already
  // built TChain with ESD trees
  InitIndex();
  InitAlignObjs();
}


//______________________________________________________________________________
AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
  fPointsFilename("AliTrackPoints.root"),
  fPointsFile(0),
  fPointsTree(0),
  fLastIndex(0),
  fArrayIndex(0),
  fIsIndexBuilt(kFALSE),
  fTrackFitter(0),
  fMinimizer(0)
{
  // Constructor in the case
  // the user provides a single ESD file
  // or a directory containing ESD files
  fESDChain = new TChain(esdtreename);
  fESDChain->Add(esdfilename);

  InitIndex();
  InitAlignObjs();
}

//______________________________________________________________________________
AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
  TObject(alignment)
{
  // Copy constructor
  // not implemented
  AliWarning("Copy constructor not implemented!");
}

//______________________________________________________________________________
AliAlignmentTracks& AliAlignmentTracks::operator= (const AliAlignmentTracks& alignment)
{
  // Asignment operator
  // not implemented
  if(this==&alignment) return *this;

  AliWarning("Asignment operator not implemented!");

  ((TObject *)this)->operator=(alignment);

  return *this;
}

//______________________________________________________________________________
AliAlignmentTracks::~AliAlignmentTracks()
{
  // Destructor
  if (fESDChain) delete fESDChain;

  DeleteIndex();
  DeleteAlignObjs();

  delete fTrackFitter;
  delete fMinimizer;

  if (fPointsFile) fPointsFile->Close();
}

//______________________________________________________________________________
void AliAlignmentTracks::AddESD(TChain *esdchain)
{
  // Add a chain with ESD files
  if (fESDChain)
    fESDChain->Add(esdchain);
  else
    fESDChain = esdchain;
}

//______________________________________________________________________________
void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
{
  // Add a single file or
  // a directory to the chain
  // with the ESD files
  if (fESDChain)
    fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
  else {
    fESDChain = new TChain(esdtreename);
    fESDChain->Add(esdfilename);
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::ProcessESD()
{
  // Analyzes and filters ESD tracks
  // Stores the selected track space points
  // into the output file

  if (!fESDChain) return;

  AliESD *esd = 0;
  fESDChain->SetBranchAddress("ESD",&esd);

  // Open the output file
  if (fPointsFilename.Data() == "") {
    AliWarning("Incorrect output filename!");
    return;
  }

  TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
  if (!pointsFile || !pointsFile->IsOpen()) {
    AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
    return;
  }

  TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
  AliTrackPointArray *array = 0;
  pointsTree->Branch("SP","AliTrackPointArray", &array);

  Int_t ievent = 0;
  while (fESDChain->GetEntry(ievent++)) {
    if (!esd) break;
    Int_t ntracks = esd->GetNumberOfTracks();
    for (Int_t itrack=0; itrack < ntracks; itrack++) {
      AliESDtrack * track = esd->GetTrack(itrack);
      if (!track) continue;
 
      UInt_t status = AliESDtrack::kITSpid; 
      status|=AliESDtrack::kTPCpid; 
      status|=AliESDtrack::kTRDpid; 
      if ((track->GetStatus() & status) != status) continue;

      if (track->GetP() < 0.5) continue;

      array = track->GetTrackPointArray();
      pointsTree->Fill();
    }
  }

  if (!pointsTree->Write()) {
    AliWarning("Can't write the tree with track point arrays!");
    return;
  }

  pointsFile->Close();
}

//______________________________________________________________________________
void AliAlignmentTracks::ProcessESD(TSelector *selector)
{
  AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
}

//______________________________________________________________________________
void AliAlignmentTracks::BuildIndex()
{
  // Build index of points tree entries
  // Used for access based on the volume IDs
  if (fIsIndexBuilt)
    ResetIndex();
  else
    fIsIndexBuilt = kTRUE;

  TFile *fPointsFile = TFile::Open(fPointsFilename);
  if (!fPointsFile || !fPointsFile->IsOpen()) {
    AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
    return;
  }
  
  //  AliTrackPointArray* array = new AliTrackPointArray;
  AliTrackPointArray* array = 0;
  TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
  if (!pointsTree) {
    AliWarning("No pointsTree found!");
    return;
  }
  pointsTree->SetBranchAddress("SP", &array);

  Int_t nArrays = pointsTree->GetEntries();
  for (Int_t iArray = 0; iArray < nArrays; iArray++)
    {
      pointsTree->GetEvent(iArray);
      if (!array) continue;
      for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
	UShort_t volId = array->GetVolumeID()[ipoint];
	Int_t modId;
	Int_t layerId = AliAlignObj::VolUIDToLayer(volId,modId)
	              - AliAlignObj::kFirstLayer;
	if (!fArrayIndex[layerId][modId]) {
	  //first entry for this volume
	  fArrayIndex[layerId][modId] = new TArrayI(1000);
	}
	else {
	  Int_t size = fArrayIndex[layerId][modId]->GetSize();
	  // If needed allocate new size
	  if (fLastIndex[layerId][modId] >= size)
	    fArrayIndex[layerId][modId]->Set(size + 1000);
	}

	// Check if the index is already filled
	Bool_t fillIndex = kTRUE;
	if (fLastIndex[layerId][modId] != 0) {
	  if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
	    fillIndex = kFALSE;
	}
	// Fill the index array and store last filled index
	if (fillIndex) {
	  (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
	  fLastIndex[layerId][modId]++;
	}
      }
    }

}

// //______________________________________________________________________________
// void AliAlignmentTracks::BuildIndexLayer(AliAlignObj::ELayerID layer)
// {
// }

// //______________________________________________________________________________
// void AliAlignmentTracks::BuildIndexVolume(UShort_t volid)
// {
// }

//______________________________________________________________________________
void AliAlignmentTracks::InitIndex()
{
  // Initialize the index arrays
  Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
  fLastIndex = new Int_t*[nLayers];
  fArrayIndex = new TArrayI**[nLayers];
  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
    fLastIndex[iLayer] = new Int_t[AliAlignObj::LayerSize(iLayer)];
    fArrayIndex[iLayer] = new TArrayI*[AliAlignObj::LayerSize(iLayer)];
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
      fLastIndex[iLayer][iModule] = 0;
      fArrayIndex[iLayer][iModule] = 0;
    }
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::ResetIndex()
{
  // Reset the value of the last filled index
  // Do not realocate memory
  for (Int_t iLayer = AliAlignObj::kFirstLayer; iLayer < AliAlignObj::kLastLayer; iLayer++) {
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
      fLastIndex[iLayer][iModule] = 0;
    }
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::DeleteIndex()
{
  // Delete the index arrays
  // Called by the destructor
  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
      if (fArrayIndex[iLayer][iModule]) {
	delete fArrayIndex[iLayer][iModule];
	fArrayIndex[iLayer][iModule] = 0;
      }
    }
    delete [] fLastIndex[iLayer];
    delete [] fArrayIndex[iLayer];
  }
  delete [] fLastIndex;
  delete [] fArrayIndex;
}

//______________________________________________________________________________
Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
{
  // Read alignment object from a file
  // To be replaced by a call to CDB
  AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));

  return kFALSE;
}

//______________________________________________________________________________
void AliAlignmentTracks::InitAlignObjs()
{
  // Initialize the alignment objects array
  Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
  fAlignObjs = new AliAlignObj**[nLayers];
  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
    fAlignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
      fAlignObjs[iLayer][iModule] = new AliAlignObjAngles;
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::ResetAlignObjs()
{
  // Reset the alignment objects array
  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
      fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::DeleteAlignObjs()
{
  // Delete the alignment objects array
  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
      if (fAlignObjs[iLayer][iModule])
	delete fAlignObjs[iLayer][iModule];
    delete [] fAlignObjs[iLayer];
  }
  delete [] fAlignObjs;
}

//______________________________________________________________________________
void AliAlignmentTracks::Align(Int_t iterations)
{
  // This method is just an example
  // how one can user AlignLayer and
  // AlignVolume methods to construct
  // a custom alignment procedure
  while (iterations > 0) {
    // First inward pass
    AlignLayer(AliAlignObj::kTPC1);
    AlignLayer(AliAlignObj::kSSD2);
    AlignLayer(AliAlignObj::kSSD1);
    AlignLayer(AliAlignObj::kSDD2);
    AlignLayer(AliAlignObj::kSDD1);
    AlignLayer(AliAlignObj::kSPD2);
    AlignLayer(AliAlignObj::kSPD1);
    // Outward pass
    AlignLayer(AliAlignObj::kSPD2);
    AlignLayer(AliAlignObj::kSDD1);
    AlignLayer(AliAlignObj::kSDD2);
    AlignLayer(AliAlignObj::kSSD1);
    AlignLayer(AliAlignObj::kSSD2);
    AlignLayer(AliAlignObj::kTPC1);
    AlignLayer(AliAlignObj::kTPC2);
    AlignLayer(AliAlignObj::kTRD1);
    AlignLayer(AliAlignObj::kTRD2);
    AlignLayer(AliAlignObj::kTRD3);
    AlignLayer(AliAlignObj::kTRD4);
    AlignLayer(AliAlignObj::kTRD5);
    AlignLayer(AliAlignObj::kTRD6);
    AlignLayer(AliAlignObj::kTOF);
    // Again inward
    AlignLayer(AliAlignObj::kTRD6);
    AlignLayer(AliAlignObj::kTRD5);
    AlignLayer(AliAlignObj::kTRD4);
    AlignLayer(AliAlignObj::kTRD3);
    AlignLayer(AliAlignObj::kTRD2);
    AlignLayer(AliAlignObj::kTRD1);
    AlignLayer(AliAlignObj::kTPC2);
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
				    AliAlignObj::ELayerID layerRangeMin,
				    AliAlignObj::ELayerID layerRangeMax,
				    Int_t iterations)
{
  // Align detector volumes within
  // a given layer.
  // Tracks are fitted only within
  // the range defined by the user.
  while (iterations > 0) {
    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
      UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
      AlignVolume(volId,layerRangeMin,layerRangeMax);
    }
    iterations--;
  }
}

//______________________________________________________________________________
void AliAlignmentTracks::AlignVolume(UShort_t volid,
				     AliAlignObj::ELayerID layerRangeMin,
				     AliAlignObj::ELayerID layerRangeMax)
{
  // Align a single detector volume.
  // Tracks are fitted only within
  // the range defined by the user.

  // First take the alignment object to be updated
  Int_t iModule;
  AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
  AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];

  // Then load only the tracks with at least one
  // space point in the volume (volid)
  BuildIndex();
  AliTrackPointArray **points;
  Int_t nArrays = LoadPoints(volid, points);
  if (nArrays == 0) return;

  AliTrackResiduals *minimizer = CreateMinimizer();
  minimizer->SetNTracks(nArrays);
  minimizer->SetAlignObj(alignObj);
  AliTrackFitter *fitter = CreateFitter();
  for (Int_t iArray = 0; iArray < nArrays; iArray++) {
    fitter->SetTrackPointArray(points[iArray], kFALSE);
    AliTrackPointArray *pVolId = 0, *pTrack = 0;
    fitter->Fit(volid,pVolId,pTrack,layerRangeMin,layerRangeMax);
    minimizer->AddTrackPointArrays(pVolId,pTrack);
  }
  minimizer->Minimize();

  UnloadPoints(nArrays, points);
}
  
//______________________________________________________________________________
Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &points)
{
  // Load track point arrays with at least
  // one space point in a given detector
  // volume (volid).
  // Use the already created tree index for
  // fast access.
  Int_t iModule;
  AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
  Int_t nArrays = fLastIndex[iLayer][iModule];

  // In case of empty index
  if (nArrays == 0) {
    points = 0;
    return 0;
  }

  if (!fPointsTree) {
    AliWarning("Tree with the space point arrays not initialized!");
    points = 0;
    return 0;
  }

  AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
  TGeoHMatrix m;
  alignObj->GetMatrix(m);
  Double_t *rot = m.GetRotationMatrix();
  Double_t *tr  = m.GetTranslation();

  AliTrackPointArray* array = 0;
  fPointsTree->SetBranchAddress("SP", &array);

  points = new AliTrackPointArray*[nArrays];
  TArrayI *index = fArrayIndex[iLayer][iModule];
  AliTrackPoint p;
  Float_t xyz[3],cov[6];
  Float_t newxyz[3];
  for (Int_t iArray = 0; iArray < nArrays; iArray++) {
    fPointsTree->GetEvent((*index)[iArray]);
    if (!array) {
      AliWarning("Wrong space point array index!");
      continue;
    }
    Int_t nPoints = array->GetNPoints();
    points[iArray] = new AliTrackPointArray(nPoints);
    for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
      array->GetPoint(p,iPoint);
      p.GetXYZ(xyz,cov);
      for (Int_t i = 0; i < 3; i++)
	newxyz[i] = tr[i] 
	          + xyz[0]*rot[3*i]
                  + xyz[1]*rot[3*i+1]
	          + xyz[2]*rot[3*i+2];
      p.SetXYZ(newxyz,cov);
      points[iArray]->AddPoint(iPoint,&p);
    }
  }

  return nArrays;
}

//______________________________________________________________________________
void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
{
  // Unload track point arrays for a given
  // detector volume
  for (Int_t iArray = 0; iArray < n; iArray++)
    delete points[iArray];
  delete [] points;
}

//______________________________________________________________________________
AliTrackFitter *AliAlignmentTracks::CreateFitter()
{
  // Check if the user has already supplied
  // a track fitter object.
  // If not, create a default one.
  if (!fTrackFitter)
    fTrackFitter = new AliTrackFitterRieman;

  return fTrackFitter;
}

//______________________________________________________________________________
AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
{
  // Check if the user has already supplied
  // a track residuals minimizer object.
  // If not, create a default one.
  if (!fMinimizer)
    fMinimizer = new AliTrackResidualsChi2;

  return fMinimizer;
}
Commit	Line	Data
98937d93	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
	27	#include "AliAlignmentTracks.h"
	28	#include "AliTrackPointArray.h"
	29	#include "AliAlignObjAngles.h"
	30	#include "AliTrackFitterRieman.h"
	31	#include "AliTrackResidualsChi2.h"
	32	#include "AliESD.h"
	33	#include "AliLog.h"
	34
	35	ClassImp(AliAlignmentTracks)
	36
	37	//______________________________________________________________________________
	38	AliAlignmentTracks::AliAlignmentTracks():
	39	fESDChain(0),
	40	fPointsFilename("AliTrackPoints.root"),
	41	fPointsFile(0),
	42	fPointsTree(0),
	43	fLastIndex(0),
	44	fArrayIndex(0),
	45	fIsIndexBuilt(kFALSE),
	46	fTrackFitter(0),
	47	fMinimizer(0)
	48	{
	49	// Default constructor
	50	InitIndex();
	51	InitAlignObjs();
	52	}
	53
	54	//______________________________________________________________________________
	55	AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
	56	fESDChain(esdchain),
	57	fPointsFilename("AliTrackPoints.root"),
	58	fPointsFile(0),
	59	fPointsTree(0),
	60	fLastIndex(0),
	61	fArrayIndex(0),
	62	fIsIndexBuilt(kFALSE),
	63	fTrackFitter(0),
	64	fMinimizer(0)
65	{
66	// Constructor in the case
67	// the user provides an already
68	// built TChain with ESD trees
69	InitIndex();
70	InitAlignObjs();
71	}
72
73
74	//______________________________________________________________________________
75	AliAlignmentTracks::AliAlignmentTracks(const char esdfilename, const char esdtreename):
76	fPointsFilename("AliTrackPoints.root"),
77	fPointsFile(0),
78	fPointsTree(0),
79	fLastIndex(0),
80	fArrayIndex(0),
81	fIsIndexBuilt(kFALSE),
82	fTrackFitter(0),
83	fMinimizer(0)
84	{
85	// Constructor in the case
86	// the user provides a single ESD file
87	// or a directory containing ESD files
88	fESDChain = new TChain(esdtreename);
89	fESDChain->Add(esdfilename);
90
91	InitIndex();
92	InitAlignObjs();
93	}
94
95	//______________________________________________________________________________
96	AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
97	TObject(alignment)
98	{
99	// Copy constructor
100	// not implemented
101	AliWarning("Copy constructor not implemented!");
102	}
103
104	//______________________________________________________________________________
105	AliAlignmentTracks& AliAlignmentTracks::operator= (const AliAlignmentTracks& alignment)
106	{
107	// Asignment operator
108	// not implemented
109	if(this==&alignment) return *this;
110
111	AliWarning("Asignment operator not implemented!");
112
113	((TObject *)this)->operator=(alignment);
114
115	return *this;
116	}
117
118	//______________________________________________________________________________
119	AliAlignmentTracks::~AliAlignmentTracks()
120	{
121	// Destructor
122	if (fESDChain) delete fESDChain;
123
124	DeleteIndex();
125	DeleteAlignObjs();
126
127	delete fTrackFitter;
128	delete fMinimizer;
129
130	if (fPointsFile) fPointsFile->Close();
131	}
132
133	//______________________________________________________________________________
134	void AliAlignmentTracks::AddESD(TChain *esdchain)
135	{
136	// Add a chain with ESD files
137	if (fESDChain)
138	fESDChain->Add(esdchain);
139	else
140	fESDChain = esdchain;
141	}
142
143	//______________________________________________________________________________
144	void AliAlignmentTracks::AddESD(const char esdfilename, const char esdtreename)
145	{
146	// Add a single file or
147	// a directory to the chain
148	// with the ESD files
149	if (fESDChain)
150	fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
151	else {
152	fESDChain = new TChain(esdtreename);
153	fESDChain->Add(esdfilename);
154	}
155	}
156
157	//______________________________________________________________________________
158	void AliAlignmentTracks::ProcessESD()
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	AliESD *esd = 0;
167	fESDChain->SetBranchAddress("ESD",&esd);
168
169	// Open the output file
170	if (fPointsFilename.Data() == "") {
171	AliWarning("Incorrect output filename!");
172	return;
173	}
174
175	TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
176	if (!pointsFile \|\| !pointsFile->IsOpen()) {
177	AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
178	return;
179	}
180
181	TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
182	AliTrackPointArray *array = 0;
183	pointsTree->Branch("SP","AliTrackPointArray", &array);
184
185	Int_t ievent = 0;
186	while (fESDChain->GetEntry(ievent++)) {
187	if (!esd) break;
188	Int_t ntracks = esd->GetNumberOfTracks();
189	for (Int_t itrack=0; itrack < ntracks; itrack++) {
190	AliESDtrack * track = esd->GetTrack(itrack);
191	if (!track) continue;
192
193	UInt_t status = AliESDtrack::kITSpid;
194	status\|=AliESDtrack::kTPCpid;
195	status\|=AliESDtrack::kTRDpid;
196	if ((track->GetStatus() & status) != status) continue;
197
198	if (track->GetP() < 0.5) continue;
199
200	array = track->GetTrackPointArray();
201	pointsTree->Fill();
202	}
203	}
204
205	if (!pointsTree->Write()) {
206	AliWarning("Can't write the tree with track point arrays!");
207	return;
208	}
209
210	pointsFile->Close();
211	}
212
213	//______________________________________________________________________________
214	void AliAlignmentTracks::ProcessESD(TSelector *selector)
215	{
216	AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
217	}
218
219	//______________________________________________________________________________
220	void AliAlignmentTracks::BuildIndex()
221	{
222	// Build index of points tree entries
223	// Used for access based on the volume IDs
224	if (fIsIndexBuilt)
225	ResetIndex();
226	else
227	fIsIndexBuilt = kTRUE;
228
229	TFile *fPointsFile = TFile::Open(fPointsFilename);
230	if (!fPointsFile \|\| !fPointsFile->IsOpen()) {
231	AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
232	return;
233	}
234
235	// AliTrackPointArray* array = new AliTrackPointArray;
236	AliTrackPointArray* array = 0;
237	TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
238	if (!pointsTree) {
239	AliWarning("No pointsTree found!");
240	return;
241	}
242	pointsTree->SetBranchAddress("SP", &array);
243
244	Int_t nArrays = pointsTree->GetEntries();
245	for (Int_t iArray = 0; iArray < nArrays; iArray++)
246	{
247	pointsTree->GetEvent(iArray);
248	if (!array) continue;
249	for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
250	UShort_t volId = array->GetVolumeID()[ipoint];
251	Int_t modId;
252	Int_t layerId = AliAlignObj::VolUIDToLayer(volId,modId)
253	- AliAlignObj::kFirstLayer;
254	if (!fArrayIndex[layerId][modId]) {
255	//first entry for this volume
256	fArrayIndex[layerId][modId] = new TArrayI(1000);
257	}
258	else {
259	Int_t size = fArrayIndex[layerId][modId]->GetSize();
260	// If needed allocate new size
261	if (fLastIndex[layerId][modId] >= size)
262	fArrayIndex[layerId][modId]->Set(size + 1000);
263	}
264
265	// Check if the index is already filled
266	Bool_t fillIndex = kTRUE;
267	if (fLastIndex[layerId][modId] != 0) {
268	if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
269	fillIndex = kFALSE;
270	}
271	// Fill the index array and store last filled index
272	if (fillIndex) {
273	(*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
274	fLastIndex[layerId][modId]++;
275	}
276	}
277	}
278
279	}
280
281	// //______________________________________________________________________________
282	// void AliAlignmentTracks::BuildIndexLayer(AliAlignObj::ELayerID layer)
283	// {
284	// }
285
286	// //______________________________________________________________________________
287	// void AliAlignmentTracks::BuildIndexVolume(UShort_t volid)
288	// {
289	// }
290
291	//______________________________________________________________________________
292	void AliAlignmentTracks::InitIndex()
293	{
294	// Initialize the index arrays
295	Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
296	fLastIndex = new Int_t*[nLayers];
297	fArrayIndex = new TArrayI**[nLayers];
298	for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
299	fLastIndex[iLayer] = new Int_t[AliAlignObj::LayerSize(iLayer)];
300	fArrayIndex[iLayer] = new TArrayI*[AliAlignObj::LayerSize(iLayer)];
301	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
302	fLastIndex[iLayer][iModule] = 0;
303	fArrayIndex[iLayer][iModule] = 0;
304	}
305	}
306	}
307
308	//______________________________________________________________________________
309	void AliAlignmentTracks::ResetIndex()
310	{
311	// Reset the value of the last filled index
312	// Do not realocate memory
313	for (Int_t iLayer = AliAlignObj::kFirstLayer; iLayer < AliAlignObj::kLastLayer; iLayer++) {
314	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
315	fLastIndex[iLayer][iModule] = 0;
316	}
317	}
318	}
319
320	//______________________________________________________________________________
321	void AliAlignmentTracks::DeleteIndex()
322	{
323	// Delete the index arrays
324	// Called by the destructor
325	for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
326	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
327	if (fArrayIndex[iLayer][iModule]) {
328	delete fArrayIndex[iLayer][iModule];
329	fArrayIndex[iLayer][iModule] = 0;
330	}
331	}
332	delete [] fLastIndex[iLayer];
333	delete [] fArrayIndex[iLayer];
334	}
335	delete [] fLastIndex;
336	delete [] fArrayIndex;
337	}
338
339	//______________________________________________________________________________
340	Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
341	{
342	// Read alignment object from a file
343	// To be replaced by a call to CDB
344	AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));
345
346	return kFALSE;
347	}
348
349	//______________________________________________________________________________
350	void AliAlignmentTracks::InitAlignObjs()
351	{
352	// Initialize the alignment objects array
353	Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
354	fAlignObjs = new AliAlignObj**[nLayers];
355	for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
356	fAlignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
357	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
358	fAlignObjs[iLayer][iModule] = new AliAlignObjAngles;
359	}
360	}
361
362	//______________________________________________________________________________
363	void AliAlignmentTracks::ResetAlignObjs()
364	{
365	// Reset the alignment objects array
366	for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
367	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
368	fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
369	}
370	}
371
372	//______________________________________________________________________________
373	void AliAlignmentTracks::DeleteAlignObjs()
374	{
375	// Delete the alignment objects array
376	for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
377	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
378	if (fAlignObjs[iLayer][iModule])
379	delete fAlignObjs[iLayer][iModule];
380	delete [] fAlignObjs[iLayer];
381	}
382	delete [] fAlignObjs;
383	}
384
385	//______________________________________________________________________________
386	void AliAlignmentTracks::Align(Int_t iterations)
387	{
388	// This method is just an example
389	// how one can user AlignLayer and
390	// AlignVolume methods to construct
391	// a custom alignment procedure
392	while (iterations > 0) {
393	// First inward pass
394	AlignLayer(AliAlignObj::kTPC1);
395	AlignLayer(AliAlignObj::kSSD2);
396	AlignLayer(AliAlignObj::kSSD1);
397	AlignLayer(AliAlignObj::kSDD2);
398	AlignLayer(AliAlignObj::kSDD1);
399	AlignLayer(AliAlignObj::kSPD2);
400	AlignLayer(AliAlignObj::kSPD1);
401	// Outward pass
402	AlignLayer(AliAlignObj::kSPD2);
403	AlignLayer(AliAlignObj::kSDD1);
404	AlignLayer(AliAlignObj::kSDD2);
405	AlignLayer(AliAlignObj::kSSD1);
406	AlignLayer(AliAlignObj::kSSD2);
407	AlignLayer(AliAlignObj::kTPC1);
408	AlignLayer(AliAlignObj::kTPC2);
409	AlignLayer(AliAlignObj::kTRD1);
410	AlignLayer(AliAlignObj::kTRD2);
411	AlignLayer(AliAlignObj::kTRD3);
412	AlignLayer(AliAlignObj::kTRD4);
413	AlignLayer(AliAlignObj::kTRD5);
414	AlignLayer(AliAlignObj::kTRD6);
415	AlignLayer(AliAlignObj::kTOF);
416	// Again inward
417	AlignLayer(AliAlignObj::kTRD6);
418	AlignLayer(AliAlignObj::kTRD5);
419	AlignLayer(AliAlignObj::kTRD4);
420	AlignLayer(AliAlignObj::kTRD3);
421	AlignLayer(AliAlignObj::kTRD2);
422	AlignLayer(AliAlignObj::kTRD1);
423	AlignLayer(AliAlignObj::kTPC2);
424	}
425	}
426
427	//______________________________________________________________________________
428	void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
429	AliAlignObj::ELayerID layerRangeMin,
430	AliAlignObj::ELayerID layerRangeMax,
431	Int_t iterations)
432	{
433	// Align detector volumes within
434	// a given layer.
435	// Tracks are fitted only within
436	// the range defined by the user.
437	while (iterations > 0) {
438	for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
439	UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
440	AlignVolume(volId,layerRangeMin,layerRangeMax);
441	}
442	iterations--;
443	}
444	}
445
446	//______________________________________________________________________________
447	void AliAlignmentTracks::AlignVolume(UShort_t volid,
448	AliAlignObj::ELayerID layerRangeMin,
449	AliAlignObj::ELayerID layerRangeMax)
450	{
451	// Align a single detector volume.
452	// Tracks are fitted only within
453	// the range defined by the user.
454
455	// First take the alignment object to be updated
456	Int_t iModule;
457	AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
458	AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
459
460	// Then load only the tracks with at least one
461	// space point in the volume (volid)
462	BuildIndex();
463	AliTrackPointArray **points;
464	Int_t nArrays = LoadPoints(volid, points);
465	if (nArrays == 0) return;
466
467	AliTrackResiduals *minimizer = CreateMinimizer();
468	minimizer->SetNTracks(nArrays);
469	minimizer->SetAlignObj(alignObj);
470	AliTrackFitter *fitter = CreateFitter();
471	for (Int_t iArray = 0; iArray < nArrays; iArray++) {
472	fitter->SetTrackPointArray(points[iArray], kFALSE);
473	AliTrackPointArray pVolId = 0, pTrack = 0;
474	fitter->Fit(volid,pVolId,pTrack,layerRangeMin,layerRangeMax);
475	minimizer->AddTrackPointArrays(pVolId,pTrack);
476	}
477	minimizer->Minimize();
478
479	UnloadPoints(nArrays, points);
480	}
481
482	//______________________________________________________________________________
483	Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &points)
484	{
485	// Load track point arrays with at least
486	// one space point in a given detector
487	// volume (volid).
488	// Use the already created tree index for
489	// fast access.
490	Int_t iModule;
491	AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
492	Int_t nArrays = fLastIndex[iLayer][iModule];
493
494	// In case of empty index
495	if (nArrays == 0) {
496	points = 0;
497	return 0;
498	}
499
500	if (!fPointsTree) {
501	AliWarning("Tree with the space point arrays not initialized!");
502	points = 0;
503	return 0;
504	}
505
506	AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
507	TGeoHMatrix m;
508	alignObj->GetMatrix(m);
509	Double_t *rot = m.GetRotationMatrix();
510	Double_t *tr = m.GetTranslation();
511
512	AliTrackPointArray* array = 0;
513	fPointsTree->SetBranchAddress("SP", &array);
514
515	points = new AliTrackPointArray*[nArrays];
516	TArrayI *index = fArrayIndex[iLayer][iModule];
517	AliTrackPoint p;
518	Float_t xyz[3],cov[6];
519	Float_t newxyz[3];
520	for (Int_t iArray = 0; iArray < nArrays; iArray++) {
521	fPointsTree->GetEvent((*index)[iArray]);
522	if (!array) {
523	AliWarning("Wrong space point array index!");
524	continue;
525	}
526	Int_t nPoints = array->GetNPoints();
527	points[iArray] = new AliTrackPointArray(nPoints);
528	for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
529	array->GetPoint(p,iPoint);
530	p.GetXYZ(xyz,cov);
531	for (Int_t i = 0; i < 3; i++)
532	newxyz[i] = tr[i]
533	+ xyz[0]rot[3i]
534	+ xyz[1]rot[3i+1]
535	+ xyz[2]rot[3i+2];
536	p.SetXYZ(newxyz,cov);
537	points[iArray]->AddPoint(iPoint,&p);
538	}
539	}
540
541	return nArrays;
542	}
543
544	//______________________________________________________________________________
545	void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
546	{
547	// Unload track point arrays for a given
548	// detector volume
549	for (Int_t iArray = 0; iArray < n; iArray++)
550	delete points[iArray];
551	delete [] points;
552	}
553
554	//______________________________________________________________________________
555	AliTrackFitter *AliAlignmentTracks::CreateFitter()
556	{
557	// Check if the user has already supplied
558	// a track fitter object.
559	// If not, create a default one.
560	if (!fTrackFitter)
561	fTrackFitter = new AliTrackFitterRieman;
562
563	return fTrackFitter;
564	}
565
566	//______________________________________________________________________________
567	AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
568	{
569	// Check if the user has already supplied
570	// a track residuals minimizer object.
571	// If not, create a default one.
572	if (!fMinimizer)
573	fMinimizer = new AliTrackResidualsChi2;
574
575	return fMinimizer;
576	}