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