X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWG1%2FAliRelAlignerKalmanArray.cxx;h=7f933748028945ac6fda95542da926d0a57ab7de;hb=045a97a78953d7c98f75d944d3cfa9395c37c1c3;hp=1f9b88a862194429d1f41588b2b676e31d1ba0d9;hpb=7371114e22ba026b9e80c7837f0cbb8745181aae;p=u%2Fmrichter%2FAliRoot.git

diff --git a/PWG1/AliRelAlignerKalmanArray.cxx b/PWG1/AliRelAlignerKalmanArray.cxx
index 1f9b88a8621..7f933748028 100644
--- a/PWG1/AliRelAlignerKalmanArray.cxx
+++ b/PWG1/AliRelAlignerKalmanArray.cxx
@@ -24,251 +24,408 @@
 //////////////////////////////////////////////////////////////////////////////
 
 #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 <AliESDEvent.h>
+#include <AliRelAlignerKalman.h>
 #include "AliRelAlignerKalmanArray.h"
 
 ClassImp(AliRelAlignerKalmanArray)
 
 //______________________________________________________________________________
 AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
-    TNamed(),
-    fArray(new TObjArray()),
-    fSaveInterval(600), //default every 10minutes
-    fTimeMatchingTolerance(20),
-    fCurrentTimeBin(0),
-    fAligner(new AliRelAlignerKalman())
+    TNamed("alignerArray","array of aligners"),
+    fT0(0),
+    fTimebinWidth(0),
+    fSize(0),
+    fOutRejSigmaOnMerge(10.),
+    fOutRejSigmaOnSmooth(1.),
+    fAlignerTemplate(),
+    fPArray(NULL)
 {
   //ctor
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(const char* name):
-    TNamed(name, name),
-    fArray(new TObjArray()),
-    fSaveInterval(600), //default every 10 minutes
-    fTimeMatchingTolerance(60),
-    fCurrentTimeBin(0),
-    fAligner(new AliRelAlignerKalman())
+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
-  fArray->SetOwner();
-  delete fArray;
-  delete fAligner;
+  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
-  //the merge is vertical, meaning matching entries in tree are merged
+  //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 (arrayFromList==this) continue;
-
-    fArray->AddAll(arrayFromList->fArray); //put all objects in one array
-
-    //do the merge
-    fArray->Sort();
-    TObjArray* tmpArray = SortedMerge(fArray);
-    tmpArray->SetOwner(kTRUE);
-
-    TObjArray* newArray = dynamic_cast<TObjArray*>(tmpArray->Clone());
-    delete fArray; //takes care of all loaded objects
-    fArray = newArray;
+    if (fT0 != arrayFromList->fT0) continue;
+    if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
+    if (fSize != arrayFromList->fSize) continue;
 
-    fArray->AddLast(arrayFromList->fAligner); //add the endofrun aligner
+    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;
+}
 
-  //TODO: this can be done better!
-  //Add own endofrun aligner and clean up
-  fArray->AddLast(fAligner);
-  fArray->Sort();
-  //TObjArray* tmpArray = SortedMerge(fArray);
-  //tmpArray->SetOwner(kTRUE);
-  //TObjArray* newArray = dynamic_cast<TObjArray*>(tmpArray->Clone());
-  //delete fArray; //takes care of all loaded objects
-  //fArray = newArray;
-
-  return fArray->GetEntriesFast();
+//______________________________________________________________________________
+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;
 }
 
 //______________________________________________________________________________
-TObjArray* AliRelAlignerKalmanArray::SortedMerge( TObjArray* input )
+AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
 {
-  //Merges the adjacent aligners if close enough
-  //input needs to be already sorted
-
-  UInt_t timeStampIn;
-  AliRelAlignerKalman* alignerIn;
-  AliRelAlignerKalman* alignerOut = dynamic_cast<AliRelAlignerKalman*>(input->At(0));
-  TObjArray* output = new TObjArray();  //empty array
-  output->AddLast(alignerOut);        //first object in: copy of first input element
-
-  timeStampIn = alignerOut->GetTimeStamp();
-  SetCurrentTimeBin( timeStampIn );
-  TIter next(input);
-  while ( (alignerIn = dynamic_cast<AliRelAlignerKalman*>(next())) )
+  //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])
   {
-    timeStampIn = alignerIn->GetTimeStamp();
-    if ( IsInCurrentTimeBin(timeStampIn) )
-    {
-      alignerOut->Merge(alignerIn);
-    }
-    else
-    {
-      alignerOut = alignerIn;
-      output->AddLast(alignerOut);
-    }//if
-    SetCurrentTimeBin( timeStampIn );
+    fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
+    fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
   }
-  return output;
+  return fPArray[tb];
 }
 
 //______________________________________________________________________________
-void AliRelAlignerKalmanArray::SetCurrentTimeBin( UInt_t timestamp )
+AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
 {
-  //set the current timebin
-  fCurrentTimeBin = TimeBin(timestamp);
+  //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;
 }
 
 //______________________________________________________________________________
-UInt_t AliRelAlignerKalmanArray::TimeBin( UInt_t timestamp ) const
+AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
 {
-  return (timestamp+(fSaveInterval/2))/fSaveInterval*fSaveInterval; //it's all integers!
+  //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;
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalmanArray::IsInCurrentTimeBin( UInt_t timestamp ) const
+void AliRelAlignerKalmanArray::ClearContents()
 {
-  //check if timestamp is within the current timebin
-  UInt_t timeDiff = (timestamp>=fCurrentTimeBin)?timestamp-fCurrentTimeBin:
-                    fCurrentTimeBin-timestamp;
-  return (timeDiff < fTimeMatchingTolerance);
+  //clear contents of array
+  for (Int_t i=0;i<fSize;i++)
+  {
+    if (fPArray[i]) delete fPArray[i];
+  }
 }
 
-////______________________________________________________________________________
-//void AliRelAlignerKalmanArray::AddESDEvent( AliESDEvent* event )
-//{
-//  //add an AliESDEvent, take care of bookkeeping
-//  if (!fAligner) return;
-//  if (event->GetRunNumber() != fAligner->GetRunNumber())
-//  {
-//    //what to do when a new run starts
-//  }
-//
-//}
-//
-////______________________________________________________________________________
-Bool_t AliRelAlignerKalmanArray::AddCosmicEvent( AliESDEvent* event )
+//______________________________________________________________________________
+AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
 {
-  if (!fAligner->AddCosmicEvent(event)) return kFALSE;
-
-  UInt_t currentTimeStamp = event->GetTimeStamp();
-  UInt_t timeFromLastBinCentre = currentTimeStamp - fCurrentTimeBin;
-  UInt_t binCentre = TimeBin(currentTimeStamp);
-  UInt_t timeFromBinCentre = currentTimeStamp-binCentre;
-  UInt_t nIntervals = timeFromLastBinCentre/fSaveInterval;
-
-  //////////////////////////////////////////////////////////////////////////////
-  //only ONE SAVE PER TIMEBIN!!!!, as close as possible to the bin center
-  //////////////////////////////////////////////////////////////////////////////
-  if ( (nIntervals == 1) &&                      //is in next time bin passed centre
-       (timeFromBinCentre < fTimeMatchingTolerance) )     //and close to it
-  {
-    AddLast(new AliRelAlignerKalman(*fAligner));
-  }
-  else if ( (nIntervals > 2) ) //TODO: don't hardwire stuff!
-  {
-    //if missed a few windows save anyway at current bin centre
-    fAligner->SetTimeStamp(binCentre);
-    AddLast(new AliRelAlignerKalman(*fAligner));
-  }
-  //////////////////////////////////////////////////////////////////////////////
+  //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];
+}
 
-  return kTRUE;
+//______________________________________________________________________________
+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];
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalmanArray::AliRelAlignerKalmanArray( const AliRelAlignerKalmanArray& in):
-    TNamed(in.GetName(), in.GetTitle()),
-    fArray(NULL),
-    fSaveInterval(in.fSaveInterval),
-    fTimeMatchingTolerance(in.fTimeMatchingTolerance),
-    fCurrentTimeBin(in.fCurrentTimeBin),
-    fAligner(new AliRelAlignerKalman(*in.fAligner))
+AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
 {
-  //copy ctor
-  fArray = static_cast<TObjArray*>(in.Clone());
+  //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];
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
+AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
 {
-  //assignment operator
-  fArray = static_cast<TObjArray*>(in.Clone());
-  fSaveInterval = in.fSaveInterval;
-  fTimeMatchingTolerance = in.fTimeMatchingTolerance;
-  return *this;
+  //return aligner in last filled slot
+  if (fSize==0) return NULL;
+  return fPArray[fSize-1];
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalmanArray::SetSaveInterval( const UInt_t s )
+void AliRelAlignerKalmanArray::Print(Option_t* option) const
 {
-  //only set if array empty
-  if (fArray->GetEntriesFast()) return kFALSE;
-  fSaveInterval = s;
-  return kTRUE;
+  // 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();
+    }
 }
 
 //______________________________________________________________________________
-Bool_t AliRelAlignerKalmanArray::SetTimeMatchingTolerance( const UInt_t m )
+Int_t AliRelAlignerKalmanArray::GetEntries() const
 {
-  //only set if array empty
-  if (fArray->GetEntriesFast()) return kFALSE;
-  fTimeMatchingTolerance = m;
-  return kTRUE;
+  //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;
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
+void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
 {
-  //mimic TObjArray::At( Int_t i )
-  return dynamic_cast<AliRelAlignerKalman*>(fArray->At(i));
+  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();
+  }
 }
 
 //______________________________________________________________________________
-void AliRelAlignerKalmanArray::AddLast( AliRelAlignerKalman* al )
+TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
 {
-  //mimic TObjArray::AddLast( TObject* obj )
-  fArray->AddLast( al );
-  SetCurrentTimeBin(al->GetTimeStamp());
+  //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;
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
+AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
 {
-  //mimic TObjArray::Last()
-  return dynamic_cast<AliRelAlignerKalman*>(fArray->Last());
+  //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;
 }
 
 //______________________________________________________________________________
-AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const
+void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
 {
-  //mimic TObjArray::operator[](Int_t)
-  return dynamic_cast<AliRelAlignerKalman*>(fArray->At(i));
+  //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;
 }