RelVal: name GNU screen after session
[u/mrichter/AliRoot.git] / 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