[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"
#include "TList.h"
#include "TBrowser.h"

ClassImp(AliRelAlignerKalmanArray)

//______________________________________________________________________________
AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
    TNamed("alignerArray","array of aligners"),
    fT0(0),
    fTimebinWidth(0),
    fSize(0),
    fOutRejSigmaOnMerge(10.),
    fOutRejSigmaOnSmooth(1.),
    fAlignerTemplate(),
    fPArray(NULL),
    fListOfGraphs(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),
    fListOfGraphs(new TList)
{
  //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;
  fListOfGraphs->SetName("graphs");
  fListOfGraphs->SetOwner(kTRUE);
}

//______________________________________________________________________________
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),
    fListOfGraphs(new TList)
{
  //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;
  }
  fListOfGraphs->SetName("graphs");
  fListOfGraphs->SetOwner(kTRUE);
}

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

//______________________________________________________________________________
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);
    }
  }
  fListOfGraphs->Delete();
  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;
  }

  Int_t n=GetEntries();
  TVectorD vx(n);
  TVectorD vy(n);
  TVectorD vex(n);
  TVectorD vey(n);
  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->SetName(graphtitle);
  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;
}

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