bug fixed for alignment, removed alignment database access from AliPMDUtility class
[u/mrichter/AliRoot.git] / 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