PWGPP/AliRelAlignerKalmanArray.cxx

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