Updating analysis task

[u/mrichter/AliRoot.git] / PWGPP / 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(&in == this) return *this;
  TNamed::operator=(in);

  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);
}