[u/mrichter/AliRoot.git] / PWG1 / AliRelAlignerKalmanArray.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.                  *
  **************************************************************************/

///////////////////////////////////////////////////////////////////////////////
//
//     Data container for relative ITS-TPC alignment analysis
//     Holds an array of AliRelAlignerKalman objects
//     and takes care of merging when processing data in parallel
//
//     Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
//
//////////////////////////////////////////////////////////////////////////////

#include <TNamed.h>
#include <TGraphErrors.h>
#include <TAxis.h>
#include <TTree.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TCollection.h>
#include <AliESDEvent.h>
#include <AliRelAlignerKalman.h>
#include "AliRelAlignerKalmanArray.h"

ClassImp(AliRelAlignerKalmanArray)

//______________________________________________________________________________
AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
    TNamed("alignerArray","array of aligners"),
    fT0(0),
    fTimebinWidth(0),
    fSize(0),
    fOutRejSigmaOnMerge(10.),
    fOutRejSigmaOnSmooth(1.),
    fAlignerTemplate(),
    fPArray(NULL)
{
  //ctor
}

//______________________________________________________________________________
AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(Int_t t0, Int_t tend, Int_t slotwidth):
    TNamed("alignerArray","array of aligners"),
    fT0(t0),
    fTimebinWidth(slotwidth),
    fSize(0),
    fOutRejSigmaOnMerge(10.),
    fOutRejSigmaOnSmooth(1.),
    fAlignerTemplate(),
    fPArray(NULL)
{
  //ctor
  if (slotwidth==0) fSize = 1;
  else fSize = (tend-t0)/slotwidth;
  fPArray = new AliRelAlignerKalman* [fSize];
  if (fPArray) for (Int_t i=0;i<fSize;i++){fPArray[i]=NULL;}//fill with zeros
  else fSize=0;
}

//______________________________________________________________________________
AliRelAlignerKalmanArray::AliRelAlignerKalmanArray( const AliRelAlignerKalmanArray& in):
    TNamed(in.GetName(), in.GetTitle()),
    fT0(in.fT0),
    fTimebinWidth(in.fTimebinWidth),
    fSize(in.fSize),
    fOutRejSigmaOnMerge(in.fOutRejSigmaOnMerge),
    fOutRejSigmaOnSmooth(in.fOutRejSigmaOnSmooth),
    fAlignerTemplate(in.fAlignerTemplate),
    fPArray(NULL)
{
  //copy ctor
  fPArray = new AliRelAlignerKalman* [fSize];
  if (!fPArray) {fSize=0;return;} //if fail
  for (Int_t i=0;i<fSize;i++)
  {
    if (in.fPArray[i]) fPArray[i]=new AliRelAlignerKalman(*(in.fPArray[i]));
    else fPArray[i]=NULL;
  }
}

//______________________________________________________________________________
AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
{
  //dtor
  ClearContents();
  delete [] fPArray;
}

//______________________________________________________________________________
void AliRelAlignerKalmanArray::SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
{
  //setup array - clears old content
  ClearContents();
  fT0=t0;
  fTimebinWidth=slotwidth;
  Int_t tmpsize;
  if (slotwidth==0) tmpsize = 1;
  else tmpsize = (tend-t0)/slotwidth;
  if (tmpsize!=fSize)
  {
    delete [] fPArray;
    fPArray=new AliRelAlignerKalman* [tmpsize];
    if (fPArray) fSize=tmpsize;
    else fSize=0;
  }
  for (Int_t i=0;i<fSize;i++){fPArray[i]=NULL;}//fill with zeros
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAlignerTemplate()
{
  //get aligner template
  return &fAlignerTemplate;
}

//______________________________________________________________________________
Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
{
  //Merge all the arrays
  //tlihe merge is vertical, meaning matching entries in tree are merged

  AliRelAlignerKalmanArray *arrayFromList;
  if (!list) return 0;
  TIter next(list);
  while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
  {
    if (fT0 != arrayFromList->fT0) continue;
    if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
    if (fSize != arrayFromList->fSize) continue;

    for (Int_t i=0; i<GetSize(); i++)
    {
      AliRelAlignerKalman* a1 = fPArray[i];
      AliRelAlignerKalman* a2 = arrayFromList->fPArray[i];
      if (a1 && a2)
      {
        a1->SetRejectOutliers(kFALSE);
        a1->SetOutRejSigma(fOutRejSigmaOnMerge); a1->Merge(a2);
      }
      else
        if (!a1 && a2) fPArray[i] = new AliRelAlignerKalman(*a2);
    }
  }
  return 0;
}

//______________________________________________________________________________
Int_t AliRelAlignerKalmanArray::Timebin( UInt_t timestamp ) const
{
  //calculate binnumber for given timestamp
  if (fTimebinWidth==0) return 0;
  Int_t slot = (timestamp-fT0)/fTimebinWidth; //it's all integers!
  return slot;
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
{
  //get the aligner for specified timestamp
  Int_t tb = Timebin(ts);
  if (tb<0) return NULL;
  if (tb>=fSize) return NULL;
  if (!fPArray) return NULL;
  if (!fPArray[tb])
  {
    fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
    fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
  }
  return fPArray[tb];
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
{
  //get the aligner for this event and set relevant info
  AliRelAlignerKalman* a = GetAligner(event->GetTimeStamp());
  if (a) a->SetRunNumber(event->GetRunNumber());
  if (a) a->SetMagField(event->GetMagneticField());
  return a;
}

//______________________________________________________________________________
AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
{
  //assignment operator
  if (fSize!=in.fSize)
  {
    //if sizes different, delete curent and make a new one with new size
    ClearContents();
    delete [] fPArray;
    fPArray = new AliRelAlignerKalman* [in.fSize];
  }

  fOutRejSigmaOnMerge=in.fOutRejSigmaOnMerge;
  fOutRejSigmaOnSmooth=in.fOutRejSigmaOnSmooth;

  if (!fPArray) fSize=0;
  else fSize = in.fSize;
  fTimebinWidth = in.fTimebinWidth;
  fT0 = in.fT0;
  for (Int_t i=0;i<fSize;i++)
  {
    if (in.fPArray[i]) fPArray[i] = new AliRelAlignerKalman(*(in.fPArray[i]));
    else fPArray[i]=NULL;
  }
  return *this;
}

//______________________________________________________________________________
void AliRelAlignerKalmanArray::ClearContents()
{
  //clear contents of array
  for (Int_t i=0;i<fSize;i++)
  {
    if (fPArray[i]) delete fPArray[i];
  }
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
{
  //mimic TObjArray::At( Int_t i )
  if (i>=fSize) {printf("AliRelAlignerKalmanArray::At index %i out of bounds, size=%i\n",i,fSize); return NULL;}
  if (i<0) return NULL;
  return fPArray[i];
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const
{
  //mimic TObjArray::operator[](Int_t)
  if (i>=fSize) {printf("AliRelAlignerKalmanArray::operator[] index %i out of bounds, size=%i\n",i,fSize); return NULL;}
  if (i<0) return NULL;
  return fPArray[i];
}

//______________________________________________________________________________
AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
{
  //mimic TObjArray::operator[](Int_t) can be used as lvalue
  if (i>=fSize||i<0) {Error("operator[]", "index %i out of bounds, size=%i\n",i,fSize); return fPArray[0];}
  return fPArray[i];
}

//______________________________________________________________________________
AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
{
  //return aligner in last filled slot
  if (fSize==0) return NULL;
  return fPArray[fSize-1];
}

//______________________________________________________________________________
void AliRelAlignerKalmanArray::Print(Option_t* option) const
{
  // print some info
  TString optionStr(option);
  printf( "%s: t0: %i, tend: %i, width: %i, size: %i, entries: %i\n",
          GetName(),
          fT0, (fT0+(fSize)*fTimebinWidth), fTimebinWidth,
          fSize, GetEntries() );
  if (optionStr.Contains("a"))
    for (Int_t i=0; i<fSize; i++)
    {
      AliRelAlignerKalman* al = fPArray[i];
      if (!al) continue;
      al->Print();
    }
}

//______________________________________________________________________________
Int_t AliRelAlignerKalmanArray::GetEntries() const
{
  //get number of filled slots
  if (!fPArray) return 0;
  Int_t entries=0;
  for (Int_t i=0; i<fSize; i++)
  {
    if (fPArray[i]) entries++;
  }
  return entries;
}

//______________________________________________________________________________
void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
{
  AliRelAlignerKalman* al = NULL;
  tree->Branch("aligner","AliRelAlignerKalman",&al);
  //fill a tree with filled slots
  for (Int_t i=0; i<fSize; i++)
  {
    al = fPArray[i];
    if (al) tree->Fill();
  }
}

//______________________________________________________________________________
TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
{
  //return a graph for the iparam-th parameter
  if (iparam>8 || iparam<0)
  {
    printf("wrong parameter number. must be from 0-8");
    return NULL;
  }

  TVectorD vx(GetEntries());
  TVectorD vy(GetEntries());
  TVectorD vex(GetEntries());
  TVectorD vey(GetEntries());
  Int_t entry=0;
  for (Int_t i=0; i<fSize; i++)
  {
    AliRelAlignerKalman* al = fPArray[i];
    if (!al) continue;
    vx(entry) = al->GetTimeStamp();
    vy(entry) = al->GetStateArr()[iparam];
    TMatrixDSym* cm = al->GetStateCov();
    vey(entry) = TMath::Sqrt((*cm)(iparam,iparam));
    entry++;
  }

  TGraphErrors* graph = new TGraphErrors(vx,vy,vex,vey);
  TString graphtitle;
  TString graphtitley;
  switch (iparam)
  {
  case 0:
    graphtitle="rotation \\psi";
    graphtitley="\\psi [rad]";
    break;
  case 1:
    graphtitle="rotation \\theta";
    graphtitley="\\theta [rad]";
    break;
  case 2:
    graphtitle="rotation \\phi";
    graphtitley="\\phi [rad]";
    break;
  case 3:
    graphtitle="shift x";
    graphtitley="x [cm]";
    break;
  case 4:
    graphtitle="shift y";
    graphtitley="y [cm]";
    break;
  case 5:
    graphtitle="shift z";
    graphtitley="z [cm]";
    break;
  case 6:
    graphtitle="TPC vd correction";
    graphtitley="correction factor []";
    break;
  case 7:
    graphtitle="TPC t0 correction";
    graphtitley="t0 correction [\\micros]";
    break;
  case 8:
    graphtitle="TPC dv/dy";
    graphtitley="dv/dy [cm/\\micros/m]";
    break;
  }
  graph->SetTitle(graphtitle);
  TAxis* xas = graph->GetXaxis();
  TAxis* yas = graph->GetYaxis();
  xas->SetTitle("time");
  xas->SetTimeDisplay(1);
  yas->SetTitle(graphtitley);
  return graph;
}

//______________________________________________________________________________
AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
{
  //return a smoothed version of the data
  AliRelAlignerKalmanArray* outputarr = new AliRelAlignerKalmanArray(fT0,
                                        fT0+fSize*fTimebinWidth,fTimebinWidth);

  AliRelAlignerKalman* tmpaligner = outputarr->GetAlignerTemplate();
  tmpaligner->SetOutRejSigma(fOutRejSigmaOnSmooth);
  //copy the first filled slot
  Int_t n=0;
  while (!fPArray[n]) {n++;}
  if (n==fSize) return NULL;
  *tmpaligner   = *fPArray[n];
  if (fPArray[n]->GetNUpdates()>10)
  (*outputarr)[n] = new AliRelAlignerKalman(*(fPArray[n]));
  //for the rest of slots use merge
  for (Int_t i=n+1;i<fSize;i++)
  {
    if (!fPArray[i]) continue;
    PropagateToTime(tmpaligner, fPArray[i]->GetTimeStamp());
    if (tmpaligner->Merge(fPArray[i]))
    (*outputarr)[i] = new AliRelAlignerKalman(*tmpaligner);
    else
    (*outputarr)[i] = new AliRelAlignerKalman(*(fPArray[i]));
  }
  
  tmpaligner->Reset(); //clean up

  return outputarr;
}

//______________________________________________________________________________
void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
{
  //propagate an aligner in time
  TMatrixDSym* cov = al->GetStateCov();
  TMatrixDSym corr(TMatrixDSym::kZero,*cov);
  UInt_t dt = (timestamp>al->GetTimeStamp())?
    timestamp-al->GetTimeStamp():al->GetTimeStamp()-timestamp;
  //the propagation matrix
  //to be added to the covariance matrix of kalman filter
  corr(6,6) = dt*1e-8;  //vdrift
  (*cov) += corr;
}
Commit	Line	Data
4a84c20d	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	//
	18	// Data container for relative ITS-TPC alignment analysis
	19	// Holds an array of AliRelAlignerKalman objects
	20	// and takes care of merging when processing data in parallel
	21	//
	22	// Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
	23	//
	24	//////////////////////////////////////////////////////////////////////////////
	25
	26	#include <TNamed.h>
51cfc50d	27	#include <TGraphErrors.h>
	28	#include <TAxis.h>
	29	#include <TTree.h>
4a84c20d	30	#include <TMath.h>
	31	#include <TObjArray.h>
	32	#include <TCollection.h>
51cfc50d	33	#include <AliESDEvent.h>
8478bf3d	34	#include <AliRelAlignerKalman.h>
4a84c20d	35	#include "AliRelAlignerKalmanArray.h"
	36
	37	ClassImp(AliRelAlignerKalmanArray)
	38
	39	//______________________________________________________________________________
	40	AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
51cfc50d	41	TNamed("alignerArray","array of aligners"),
	42	fT0(0),
	43	fTimebinWidth(0),
	44	fSize(0),
	45	fOutRejSigmaOnMerge(10.),
	46	fOutRejSigmaOnSmooth(1.),
8478bf3d	47	fAlignerTemplate(),
51cfc50d	48	fPArray(NULL)
4a84c20d	49	{
	50	//ctor
	51	}
	52
	53	//______________________________________________________________________________
51cfc50d	54	AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(Int_t t0, Int_t tend, Int_t slotwidth):
	55	TNamed("alignerArray","array of aligners"),
	56	fT0(t0),
	57	fTimebinWidth(slotwidth),
	58	fSize(0),
	59	fOutRejSigmaOnMerge(10.),
	60	fOutRejSigmaOnSmooth(1.),
8478bf3d	61	fAlignerTemplate(),
51cfc50d	62	fPArray(NULL)
4a84c20d	63	{
4a84c20d	64	//ctor
51cfc50d	65	if (slotwidth==0) fSize = 1;
	66	else fSize = (tend-t0)/slotwidth;
	67	fPArray = new AliRelAlignerKalman* [fSize];
	68	if (fPArray) for (Int_t i=0;i<fSize;i++){fPArray[i]=NULL;}//fill with zeros
	69	else fSize=0;
	70	}
	71
	72	//______________________________________________________________________________
	73	AliRelAlignerKalmanArray::AliRelAlignerKalmanArray( const AliRelAlignerKalmanArray& in):
	74	TNamed(in.GetName(), in.GetTitle()),
	75	fT0(in.fT0),
	76	fTimebinWidth(in.fTimebinWidth),
	77	fSize(in.fSize),
	78	fOutRejSigmaOnMerge(in.fOutRejSigmaOnMerge),
	79	fOutRejSigmaOnSmooth(in.fOutRejSigmaOnSmooth),
8478bf3d	80	fAlignerTemplate(in.fAlignerTemplate),
51cfc50d	81	fPArray(NULL)
	82	{
	83	//copy ctor
51cfc50d	84	fPArray = new AliRelAlignerKalman* [fSize];
	85	if (!fPArray) {fSize=0;return;} //if fail
	86	for (Int_t i=0;i<fSize;i++)
	87	{
	88	if (in.fPArray[i]) fPArray[i]=new AliRelAlignerKalman(*(in.fPArray[i]));
	89	else fPArray[i]=NULL;
	90	}
4a84c20d	91	}
	92
	93	//______________________________________________________________________________
	94	AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
	95	{
	96	//dtor
51cfc50d	97	ClearContents();
	98	delete [] fPArray;
	99	}
	100
	101	//______________________________________________________________________________
	102	void AliRelAlignerKalmanArray::SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
	103	{
	104	//setup array - clears old content
	105	ClearContents();
	106	fT0=t0;
	107	fTimebinWidth=slotwidth;
	108	Int_t tmpsize;
	109	if (slotwidth==0) tmpsize = 1;
	110	else tmpsize = (tend-t0)/slotwidth;
	111	if (tmpsize!=fSize)
	112	{
	113	delete [] fPArray;
	114	fPArray=new AliRelAlignerKalman* [tmpsize];
	115	if (fPArray) fSize=tmpsize;
	116	else fSize=0;
	117	}
	118	for (Int_t i=0;i<fSize;i++){fPArray[i]=NULL;}//fill with zeros
	119	}
	120
	121	//______________________________________________________________________________
	122	AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAlignerTemplate()
	123	{
8478bf3d	124	//get aligner template
8478bf3d	125	return &fAlignerTemplate;
4a84c20d	126	}
	127
	128	//______________________________________________________________________________
	129	Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
	130	{
	131	//Merge all the arrays
51cfc50d	132	//tlihe merge is vertical, meaning matching entries in tree are merged
4a84c20d	133
	134	AliRelAlignerKalmanArray *arrayFromList;
	135	if (!list) return 0;
	136	TIter next(list);
7371114e	137	while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
4a84c20d	138	{
51cfc50d	139	if (fT0 != arrayFromList->fT0) continue;
	140	if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
	141	if (fSize != arrayFromList->fSize) continue;
4a84c20d	142
51cfc50d	143	for (Int_t i=0; i<GetSize(); i++)
	144	{
	145	AliRelAlignerKalman* a1 = fPArray[i];
	146	AliRelAlignerKalman* a2 = arrayFromList->fPArray[i];
	147	if (a1 && a2)
	148	{
	149	a1->SetRejectOutliers(kFALSE);
	150	a1->SetOutRejSigma(fOutRejSigmaOnMerge); a1->Merge(a2);
	151	}
	152	else
	153	if (!a1 && a2) fPArray[i] = new AliRelAlignerKalman(*a2);
	154	}
4a84c20d	155	}
51cfc50d	156	return 0;
51cfc50d	157	}
4a84c20d	158
51cfc50d	159	//______________________________________________________________________________
	160	Int_t AliRelAlignerKalmanArray::Timebin( UInt_t timestamp ) const
	161	{
	162	//calculate binnumber for given timestamp
	163	if (fTimebinWidth==0) return 0;
	164	Int_t slot = (timestamp-fT0)/fTimebinWidth; //it's all integers!
	165	return slot;
4a84c20d	166	}
	167
	168	//______________________________________________________________________________
51cfc50d	169	AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
4a84c20d	170	{
51cfc50d	171	//get the aligner for specified timestamp
	172	Int_t tb = Timebin(ts);
	173	if (tb<0) return NULL;
	174	if (tb>=fSize) return NULL;
	175	if (!fPArray) return NULL;
	176	if (!fPArray[tb])
4a84c20d	177	{
8478bf3d	178	fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
51cfc50d	179	fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
4a84c20d	180	}
51cfc50d	181	return fPArray[tb];
4a84c20d	182	}
	183
	184	//______________________________________________________________________________
51cfc50d	185	AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
4a84c20d	186	{
51cfc50d	187	//get the aligner for this event and set relevant info
	188	AliRelAlignerKalman* a = GetAligner(event->GetTimeStamp());
	189	if (a) a->SetRunNumber(event->GetRunNumber());
	190	if (a) a->SetMagField(event->GetMagneticField());
	191	return a;
4a84c20d	192	}
	193
	194	//______________________________________________________________________________
51cfc50d	195	AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
4a84c20d	196	{
51cfc50d	197	//assignment operator
	198	if (fSize!=in.fSize)
	199	{
	200	//if sizes different, delete curent and make a new one with new size
	201	ClearContents();
	202	delete [] fPArray;
	203	fPArray = new AliRelAlignerKalman* [in.fSize];
	204	}
	205
	206	fOutRejSigmaOnMerge=in.fOutRejSigmaOnMerge;
	207	fOutRejSigmaOnSmooth=in.fOutRejSigmaOnSmooth;
	208
	209	if (!fPArray) fSize=0;
	210	else fSize = in.fSize;
	211	fTimebinWidth = in.fTimebinWidth;
	212	fT0 = in.fT0;
	213	for (Int_t i=0;i<fSize;i++)
	214	{
	215	if (in.fPArray[i]) fPArray[i] = new AliRelAlignerKalman(*(in.fPArray[i]));
	216	else fPArray[i]=NULL;
	217	}
	218	return *this;
4a84c20d	219	}
	220
	221	//______________________________________________________________________________
51cfc50d	222	void AliRelAlignerKalmanArray::ClearContents()
4a84c20d	223	{
51cfc50d	224	//clear contents of array
	225	for (Int_t i=0;i<fSize;i++)
	226	{
	227	if (fPArray[i]) delete fPArray[i];
	228	}
4a84c20d	229	}
4a84c20d	230
51cfc50d	231	//______________________________________________________________________________
51cfc50d	232	AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
4a84c20d	233	{
51cfc50d	234	//mimic TObjArray::At( Int_t i )
	235	if (i>=fSize) {printf("AliRelAlignerKalmanArray::At index %i out of bounds, size=%i\n",i,fSize); return NULL;}
	236	if (i<0) return NULL;
	237	return fPArray[i];
	238	}
4a84c20d	239
51cfc50d	240	//______________________________________________________________________________
	241	AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const
	242	{
	243	//mimic TObjArray::operator[](Int_t)
	244	if (i>=fSize) {printf("AliRelAlignerKalmanArray::operator[] index %i out of bounds, size=%i\n",i,fSize); return NULL;}
	245	if (i<0) return NULL;
	246	return fPArray[i];
4a84c20d	247	}
	248
	249	//______________________________________________________________________________
51cfc50d	250	AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
4a84c20d	251	{
51cfc50d	252	//mimic TObjArray::operator[](Int_t) can be used as lvalue
	253	if (i>=fSize\|\|i<0) {Error("operator[]", "index %i out of bounds, size=%i\n",i,fSize); return fPArray[0];}
	254	return fPArray[i];
4a84c20d	255	}
	256
	257	//______________________________________________________________________________
51cfc50d	258	AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
4a84c20d	259	{
51cfc50d	260	//return aligner in last filled slot
	261	if (fSize==0) return NULL;
	262	return fPArray[fSize-1];
4a84c20d	263	}
	264
	265	//______________________________________________________________________________
51cfc50d	266	void AliRelAlignerKalmanArray::Print(Option_t* option) const
4a84c20d	267	{
51cfc50d	268	// print some info
	269	TString optionStr(option);
	270	printf( "%s: t0: %i, tend: %i, width: %i, size: %i, entries: %i\n",
	271	GetName(),
	272	fT0, (fT0+(fSize)*fTimebinWidth), fTimebinWidth,
	273	fSize, GetEntries() );
	274	if (optionStr.Contains("a"))
	275	for (Int_t i=0; i<fSize; i++)
	276	{
	277	AliRelAlignerKalman* al = fPArray[i];
	278	if (!al) continue;
	279	al->Print();
	280	}
4a84c20d	281	}
	282
	283	//______________________________________________________________________________
51cfc50d	284	Int_t AliRelAlignerKalmanArray::GetEntries() const
4a84c20d	285	{
51cfc50d	286	//get number of filled slots
	287	if (!fPArray) return 0;
	288	Int_t entries=0;
	289	for (Int_t i=0; i<fSize; i++)
	290	{
	291	if (fPArray[i]) entries++;
	292	}
	293	return entries;
4a84c20d	294	}
	295
	296	//______________________________________________________________________________
51cfc50d	297	void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
4a84c20d	298	{
51cfc50d	299	AliRelAlignerKalman* al = NULL;
	300	tree->Branch("aligner","AliRelAlignerKalman",&al);
	301	//fill a tree with filled slots
	302	for (Int_t i=0; i<fSize; i++)
	303	{
	304	al = fPArray[i];
	305	if (al) tree->Fill();
	306	}
4a84c20d	307	}
	308
	309	//______________________________________________________________________________
51cfc50d	310	TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
4a84c20d	311	{
51cfc50d	312	//return a graph for the iparam-th parameter
	313	if (iparam>8 \|\| iparam<0)
	314	{
	315	printf("wrong parameter number. must be from 0-8");
	316	return NULL;
	317	}
	318
	319	TVectorD vx(GetEntries());
	320	TVectorD vy(GetEntries());
	321	TVectorD vex(GetEntries());
	322	TVectorD vey(GetEntries());
	323	Int_t entry=0;
	324	for (Int_t i=0; i<fSize; i++)
	325	{
	326	AliRelAlignerKalman* al = fPArray[i];
	327	if (!al) continue;
	328	vx(entry) = al->GetTimeStamp();
	329	vy(entry) = al->GetStateArr()[iparam];
	330	TMatrixDSym* cm = al->GetStateCov();
	331	vey(entry) = TMath::Sqrt((*cm)(iparam,iparam));
	332	entry++;
	333	}
	334
	335	TGraphErrors* graph = new TGraphErrors(vx,vy,vex,vey);
	336	TString graphtitle;
	337	TString graphtitley;
	338	switch (iparam)
	339	{
	340	case 0:
	341	graphtitle="rotation \\psi";
	342	graphtitley="\\psi [rad]";
	343	break;
	344	case 1:
	345	graphtitle="rotation \\theta";
	346	graphtitley="\\theta [rad]";
	347	break;
	348	case 2:
	349	graphtitle="rotation \\phi";
	350	graphtitley="\\phi [rad]";
	351	break;
	352	case 3:
	353	graphtitle="shift x";
	354	graphtitley="x [cm]";
	355	break;
	356	case 4:
	357	graphtitle="shift y";
	358	graphtitley="y [cm]";
	359	break;
	360	case 5:
	361	graphtitle="shift z";
	362	graphtitley="z [cm]";
	363	break;
	364	case 6:
	365	graphtitle="TPC vd correction";
	366	graphtitley="correction factor []";
	367	break;
	368	case 7:
	369	graphtitle="TPC t0 correction";
	370	graphtitley="t0 correction [\\micros]";
	371	break;
	372	case 8:
	373	graphtitle="TPC dv/dy";
	374	graphtitley="dv/dy [cm/\\micros/m]";
	375	break;
376	}
377	graph->SetTitle(graphtitle);
378	TAxis* xas = graph->GetXaxis();
379	TAxis* yas = graph->GetYaxis();
380	xas->SetTitle("time");
381	xas->SetTimeDisplay(1);
382	yas->SetTitle(graphtitley);
383	return graph;
4a84c20d	384	}
	385
	386	//______________________________________________________________________________
51cfc50d	387	AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
4a84c20d	388	{
51cfc50d	389	//return a smoothed version of the data
	390	AliRelAlignerKalmanArray* outputarr = new AliRelAlignerKalmanArray(fT0,
	391	fT0+fSize*fTimebinWidth,fTimebinWidth);
	392
	393	AliRelAlignerKalman* tmpaligner = outputarr->GetAlignerTemplate();
	394	tmpaligner->SetOutRejSigma(fOutRejSigmaOnSmooth);
	395	//copy the first filled slot
	396	Int_t n=0;
	397	while (!fPArray[n]) {n++;}
	398	if (n==fSize) return NULL;
	399	tmpaligner = fPArray[n];
	400	if (fPArray[n]->GetNUpdates()>10)
	401	(outputarr)[n] = new AliRelAlignerKalman((fPArray[n]));
	402	//for the rest of slots use merge
	403	for (Int_t i=n+1;i<fSize;i++)
	404	{
	405	if (!fPArray[i]) continue;
	406	PropagateToTime(tmpaligner, fPArray[i]->GetTimeStamp());
	407	if (tmpaligner->Merge(fPArray[i]))
	408	(outputarr)[i] = new AliRelAlignerKalman(tmpaligner);
	409	else
	410	(outputarr)[i] = new AliRelAlignerKalman((fPArray[i]));
	411	}
	412
	413	tmpaligner->Reset(); //clean up
	414
	415	return outputarr;
4a84c20d	416	}
	417
	418	//______________________________________________________________________________
51cfc50d	419	void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
4a84c20d	420	{
51cfc50d	421	//propagate an aligner in time
	422	TMatrixDSym* cov = al->GetStateCov();
	423	TMatrixDSym corr(TMatrixDSym::kZero,*cov);
	424	UInt_t dt = (timestamp>al->GetTimeStamp())?
	425	timestamp-al->GetTimeStamp():al->GetTimeStamp()-timestamp;
	426	//the propagation matrix
	427	//to be added to the covariance matrix of kalman filter
	428	corr(6,6) = dt*1e-8; //vdrift
	429	(*cov) += corr;
4a84c20d	430	}
4a84c20d	431