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