PWG1/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
  37 ClassImp(AliRelAlignerKalmanArray)
  38
  39 //______________________________________________________________________________
  40 AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
  41     TNamed("alignerArray","array of aligners"),
  42     fT0(0),
  43     fTimebinWidth(0),
  44     fSize(0),
  45     fOutRejSigmaOnMerge(10.),
  46     fOutRejSigmaOnSmooth(1.),
  47     fAlignerTemplate(),
  48     fPArray(NULL)
  49 {
  50   //ctor
  51 }
  52
  53 //______________________________________________________________________________
  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.),
  61     fAlignerTemplate(),
  62     fPArray(NULL)
  63 {
  64   //ctor
  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),
  80     fAlignerTemplate(in.fAlignerTemplate),
  81     fPArray(NULL)
  82 {
  83   //copy ctor
  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   }
  91 }
  92
  93 //______________________________________________________________________________
  94 AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
  95 {
  96   //dtor
  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 {
 124   //get aligner template
 125   return &fAlignerTemplate;
 126 }
 127
 128 //______________________________________________________________________________
 129 Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
 130 {
 131   //Merge all the arrays
 132   //tlihe merge is vertical, meaning matching entries in tree are merged
 133
 134   AliRelAlignerKalmanArray *arrayFromList;
 135   if (!list) return 0;
 136   TIter next(list);
 137   while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
 138   {
 139     if (fT0 != arrayFromList->fT0) continue;
 140     if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
 141     if (fSize != arrayFromList->fSize) continue;
 142
 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     }
 155   }
 156   return 0;
 157 }
 158
 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;
 166 }
 167
 168 //______________________________________________________________________________
 169 AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
 170 {
 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])
 177   {
 178     fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
 179     fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
 180   }
 181   return fPArray[tb];
 182 }
 183
 184 //______________________________________________________________________________
 185 AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
 186 {
 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;
 192 }
 193
 194 //______________________________________________________________________________
 195 AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
 196 {
 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;
 219 }
 220
 221 //______________________________________________________________________________
 222 void AliRelAlignerKalmanArray::ClearContents()
 223 {
 224   //clear contents of array
 225   for (Int_t i=0;i<fSize;i++)
 226   {
 227     if (fPArray[i]) delete fPArray[i];
 228   }
 229 }
 230
 231 //______________________________________________________________________________
 232 AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
 233 {
 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 }
 239
 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];
 247 }
 248
 249 //______________________________________________________________________________
 250 AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
 251 {
 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];
 255 }
 256
 257 //______________________________________________________________________________
 258 AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
 259 {
 260   //return aligner in last filled slot
 261   if (fSize==0) return NULL;
 262   return fPArray[fSize-1];
 263 }
 264
 265 //______________________________________________________________________________
 266 void AliRelAlignerKalmanArray::Print(Option_t* option) const
 267 {
 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     }
 281 }
 282
 283 //______________________________________________________________________________
 284 Int_t AliRelAlignerKalmanArray::GetEntries() const
 285 {
 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;
 294 }
 295
 296 //______________________________________________________________________________
 297 void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
 298 {
 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   }
 307 }
 308
 309 //______________________________________________________________________________
 310 TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
 311 {
 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;
 384 }
 385
 386 //______________________________________________________________________________
 387 AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
 388 {
 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;
 416 }
 417
 418 //______________________________________________________________________________
 419 void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
 420 {
 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;
 430 }
 431