AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliRelAlignerKalmanArray.cxx
Go to the documentation of this file.
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 
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 //
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>
36 #include "TList.h"
37 #include "TBrowser.h"
38 
40 
41 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
104 {
105  //dtor
106  ClearContents();
107  delete [] fPArray;
108  delete fListOfGraphs;
109 }
110 
111 //______________________________________________________________________________
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 //______________________________________________________________________________
133 {
134  //get aligner template
135  return &fAlignerTemplate;
136 }
137 
138 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 
222 
223  if (!fPArray) fSize=0;
224  else fSize = in.fSize;
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
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 //______________________________________________________________________________
404 {
405  //return a smoothed version of the data
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 //______________________________________________________________________________
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 
Int_t Timebin(UInt_t timestamp) const
AliRelAlignerKalman * operator[](Int_t i) const
AliRelAlignerKalman * GetAligner(UInt_t timestamp)
long long Long64_t
Definition: External.C:43
Long64_t Merge(TCollection *list)
AliRelAlignerKalman ** fPArray
void SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
virtual void Print(Option_t *option="") const
AliRelAlignerKalmanArray * MakeSmoothArray() const
TGraphErrors * MakeGraph(Int_t iparam) const
const Double_t kZero
Definition: trigEffQA.C:44
AliRelAlignerKalman fAlignerTemplate
void FillTree(TTree *tree) const
void PropagateToTime(AliRelAlignerKalman *al, UInt_t timestamp) const
AliRelAlignerKalman * At(Int_t i) const
AliRelAlignerKalmanArray & operator=(const AliRelAlignerKalmanArray &a)
AliRelAlignerKalman * GetAlignerTemplate()
AliRelAlignerKalman * Last() const
const char Option_t
Definition: External.C:48