fix coverity
[u/mrichter/AliRoot.git] / PWG1 / AliRelAlignerKalmanArray.cxx
CommitLineData
4a84c20d 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>
51cfc50d 27#include <TGraphErrors.h>
28#include <TAxis.h>
29#include <TTree.h>
4a84c20d 30#include <TMath.h>
31#include <TObjArray.h>
32#include <TCollection.h>
51cfc50d 33#include <AliESDEvent.h>
8478bf3d 34#include <AliRelAlignerKalman.h>
4a84c20d 35#include "AliRelAlignerKalmanArray.h"
63f4d30f 36#include "TList.h"
37#include "TBrowser.h"
4a84c20d 38
39ClassImp(AliRelAlignerKalmanArray)
40
41//______________________________________________________________________________
42AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
51cfc50d 43 TNamed("alignerArray","array of aligners"),
44 fT0(0),
45 fTimebinWidth(0),
46 fSize(0),
47 fOutRejSigmaOnMerge(10.),
48 fOutRejSigmaOnSmooth(1.),
8478bf3d 49 fAlignerTemplate(),
63f4d30f 50 fPArray(NULL),
51 fListOfGraphs(NULL)
4a84c20d 52{
53 //ctor
54}
55
56//______________________________________________________________________________
51cfc50d 57AliRelAlignerKalmanArray::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.),
8478bf3d 64 fAlignerTemplate(),
63f4d30f 65 fPArray(NULL),
66 fListOfGraphs(new TList)
4a84c20d 67{
68 //ctor
51cfc50d 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;
63f4d30f 74 fListOfGraphs->SetName("graphs");
75 fListOfGraphs->SetOwner(kTRUE);
51cfc50d 76}
77
78//______________________________________________________________________________
79AliRelAlignerKalmanArray::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),
8478bf3d 86 fAlignerTemplate(in.fAlignerTemplate),
63f4d30f 87 fPArray(NULL),
88 fListOfGraphs(new TList)
51cfc50d 89{
90 //copy ctor
51cfc50d 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 }
63f4d30f 98 fListOfGraphs->SetName("graphs");
99 fListOfGraphs->SetOwner(kTRUE);
4a84c20d 100}
101
102//______________________________________________________________________________
103AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
104{
105 //dtor
51cfc50d 106 ClearContents();
107 delete [] fPArray;
63f4d30f 108 delete fListOfGraphs;
51cfc50d 109}
110
111//______________________________________________________________________________
112void 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//______________________________________________________________________________
132AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAlignerTemplate()
133{
8478bf3d 134 //get aligner template
135 return &fAlignerTemplate;
4a84c20d 136}
137
138//______________________________________________________________________________
139Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
140{
141 //Merge all the arrays
51cfc50d 142 //tlihe merge is vertical, meaning matching entries in tree are merged
4a84c20d 143
144 AliRelAlignerKalmanArray *arrayFromList;
145 if (!list) return 0;
146 TIter next(list);
7371114e 147 while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
4a84c20d 148 {
51cfc50d 149 if (fT0 != arrayFromList->fT0) continue;
150 if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
151 if (fSize != arrayFromList->fSize) continue;
4a84c20d 152
51cfc50d 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 }
4a84c20d 165 }
63f4d30f 166 fListOfGraphs->Delete();
51cfc50d 167 return 0;
168}
4a84c20d 169
51cfc50d 170//______________________________________________________________________________
171Int_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;
4a84c20d 177}
178
179//______________________________________________________________________________
51cfc50d 180AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
4a84c20d 181{
51cfc50d 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])
4a84c20d 188 {
8478bf3d 189 fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
51cfc50d 190 fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
4a84c20d 191 }
51cfc50d 192 return fPArray[tb];
4a84c20d 193}
194
195//______________________________________________________________________________
51cfc50d 196AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
4a84c20d 197{
51cfc50d 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;
4a84c20d 203}
204
205//______________________________________________________________________________
51cfc50d 206AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
4a84c20d 207{
51cfc50d 208 //assignment operator
209 if (fSize!=in.fSize)
210 {
211 //if sizes different, delete curent and make a new one with new size
212 ClearContents();
213 delete [] fPArray;
214 fPArray = new AliRelAlignerKalman* [in.fSize];
215 }
216
217 fOutRejSigmaOnMerge=in.fOutRejSigmaOnMerge;
218 fOutRejSigmaOnSmooth=in.fOutRejSigmaOnSmooth;
219
220 if (!fPArray) fSize=0;
221 else fSize = in.fSize;
222 fTimebinWidth = in.fTimebinWidth;
223 fT0 = in.fT0;
224 for (Int_t i=0;i<fSize;i++)
225 {
226 if (in.fPArray[i]) fPArray[i] = new AliRelAlignerKalman(*(in.fPArray[i]));
227 else fPArray[i]=NULL;
228 }
229 return *this;
4a84c20d 230}
231
232//______________________________________________________________________________
51cfc50d 233void AliRelAlignerKalmanArray::ClearContents()
4a84c20d 234{
51cfc50d 235 //clear contents of array
236 for (Int_t i=0;i<fSize;i++)
237 {
238 if (fPArray[i]) delete fPArray[i];
239 }
4a84c20d 240}
241
51cfc50d 242//______________________________________________________________________________
243AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
4a84c20d 244{
51cfc50d 245 //mimic TObjArray::At( Int_t i )
246 if (i>=fSize) {printf("AliRelAlignerKalmanArray::At index %i out of bounds, size=%i\n",i,fSize); return NULL;}
247 if (i<0) return NULL;
248 return fPArray[i];
249}
4a84c20d 250
51cfc50d 251//______________________________________________________________________________
252AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const
253{
254 //mimic TObjArray::operator[](Int_t)
255 if (i>=fSize) {printf("AliRelAlignerKalmanArray::operator[] index %i out of bounds, size=%i\n",i,fSize); return NULL;}
256 if (i<0) return NULL;
257 return fPArray[i];
4a84c20d 258}
259
260//______________________________________________________________________________
51cfc50d 261AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
4a84c20d 262{
51cfc50d 263 //mimic TObjArray::operator[](Int_t) can be used as lvalue
264 if (i>=fSize||i<0) {Error("operator[]", "index %i out of bounds, size=%i\n",i,fSize); return fPArray[0];}
265 return fPArray[i];
4a84c20d 266}
267
268//______________________________________________________________________________
51cfc50d 269AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
4a84c20d 270{
51cfc50d 271 //return aligner in last filled slot
272 if (fSize==0) return NULL;
273 return fPArray[fSize-1];
4a84c20d 274}
275
276//______________________________________________________________________________
51cfc50d 277void AliRelAlignerKalmanArray::Print(Option_t* option) const
4a84c20d 278{
51cfc50d 279 // print some info
280 TString optionStr(option);
281 printf( "%s: t0: %i, tend: %i, width: %i, size: %i, entries: %i\n",
282 GetName(),
283 fT0, (fT0+(fSize)*fTimebinWidth), fTimebinWidth,
284 fSize, GetEntries() );
285 if (optionStr.Contains("a"))
286 for (Int_t i=0; i<fSize; i++)
287 {
288 AliRelAlignerKalman* al = fPArray[i];
289 if (!al) continue;
290 al->Print();
291 }
4a84c20d 292}
293
294//______________________________________________________________________________
51cfc50d 295Int_t AliRelAlignerKalmanArray::GetEntries() const
4a84c20d 296{
51cfc50d 297 //get number of filled slots
298 if (!fPArray) return 0;
299 Int_t entries=0;
300 for (Int_t i=0; i<fSize; i++)
301 {
302 if (fPArray[i]) entries++;
303 }
304 return entries;
4a84c20d 305}
306
307//______________________________________________________________________________
51cfc50d 308void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
4a84c20d 309{
51cfc50d 310 AliRelAlignerKalman* al = NULL;
311 tree->Branch("aligner","AliRelAlignerKalman",&al);
312 //fill a tree with filled slots
313 for (Int_t i=0; i<fSize; i++)
314 {
315 al = fPArray[i];
316 if (al) tree->Fill();
317 }
4a84c20d 318}
319
320//______________________________________________________________________________
51cfc50d 321TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
4a84c20d 322{
51cfc50d 323 //return a graph for the iparam-th parameter
324 if (iparam>8 || iparam<0)
325 {
326 printf("wrong parameter number. must be from 0-8");
327 return NULL;
328 }
329
63f4d30f 330 Int_t n=GetEntries();
331 TVectorD vx(n);
332 TVectorD vy(n);
333 TVectorD vex(n);
334 TVectorD vey(n);
51cfc50d 335 Int_t entry=0;
336 for (Int_t i=0; i<fSize; i++)
337 {
338 AliRelAlignerKalman* al = fPArray[i];
339 if (!al) continue;
340 vx(entry) = al->GetTimeStamp();
341 vy(entry) = al->GetStateArr()[iparam];
342 TMatrixDSym* cm = al->GetStateCov();
343 vey(entry) = TMath::Sqrt((*cm)(iparam,iparam));
344 entry++;
345 }
346
347 TGraphErrors* graph = new TGraphErrors(vx,vy,vex,vey);
348 TString graphtitle;
349 TString graphtitley;
350 switch (iparam)
351 {
352 case 0:
353 graphtitle="rotation \\psi";
354 graphtitley="\\psi [rad]";
355 break;
356 case 1:
357 graphtitle="rotation \\theta";
358 graphtitley="\\theta [rad]";
359 break;
360 case 2:
361 graphtitle="rotation \\phi";
362 graphtitley="\\phi [rad]";
363 break;
364 case 3:
365 graphtitle="shift x";
366 graphtitley="x [cm]";
367 break;
368 case 4:
369 graphtitle="shift y";
370 graphtitley="y [cm]";
371 break;
372 case 5:
373 graphtitle="shift z";
374 graphtitley="z [cm]";
375 break;
376 case 6:
377 graphtitle="TPC vd correction";
378 graphtitley="correction factor []";
379 break;
380 case 7:
381 graphtitle="TPC t0 correction";
382 graphtitley="t0 correction [\\micros]";
383 break;
384 case 8:
385 graphtitle="TPC dv/dy";
386 graphtitley="dv/dy [cm/\\micros/m]";
387 break;
388 }
63f4d30f 389 graph->SetName(graphtitle);
51cfc50d 390 graph->SetTitle(graphtitle);
391 TAxis* xas = graph->GetXaxis();
392 TAxis* yas = graph->GetYaxis();
393 xas->SetTitle("time");
394 xas->SetTimeDisplay(1);
395 yas->SetTitle(graphtitley);
396 return graph;
4a84c20d 397}
398
399//______________________________________________________________________________
51cfc50d 400AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
4a84c20d 401{
51cfc50d 402 //return a smoothed version of the data
403 AliRelAlignerKalmanArray* outputarr = new AliRelAlignerKalmanArray(fT0,
404 fT0+fSize*fTimebinWidth,fTimebinWidth);
405
406 AliRelAlignerKalman* tmpaligner = outputarr->GetAlignerTemplate();
407 tmpaligner->SetOutRejSigma(fOutRejSigmaOnSmooth);
408 //copy the first filled slot
409 Int_t n=0;
410 while (!fPArray[n]) {n++;}
411 if (n==fSize) return NULL;
412 *tmpaligner = *fPArray[n];
413 if (fPArray[n]->GetNUpdates()>10)
414 (*outputarr)[n] = new AliRelAlignerKalman(*(fPArray[n]));
415 //for the rest of slots use merge
416 for (Int_t i=n+1;i<fSize;i++)
417 {
418 if (!fPArray[i]) continue;
419 PropagateToTime(tmpaligner, fPArray[i]->GetTimeStamp());
420 if (tmpaligner->Merge(fPArray[i]))
421 (*outputarr)[i] = new AliRelAlignerKalman(*tmpaligner);
422 else
423 (*outputarr)[i] = new AliRelAlignerKalman(*(fPArray[i]));
424 }
425
426 tmpaligner->Reset(); //clean up
427
428 return outputarr;
4a84c20d 429}
430
431//______________________________________________________________________________
51cfc50d 432void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
4a84c20d 433{
51cfc50d 434 //propagate an aligner in time
435 TMatrixDSym* cov = al->GetStateCov();
436 TMatrixDSym corr(TMatrixDSym::kZero,*cov);
437 UInt_t dt = (timestamp>al->GetTimeStamp())?
438 timestamp-al->GetTimeStamp():al->GetTimeStamp()-timestamp;
439 //the propagation matrix
440 //to be added to the covariance matrix of kalman filter
441 corr(6,6) = dt*1e-8; //vdrift
442 (*cov) += corr;
4a84c20d 443}
444
63f4d30f 445//______________________________________________________________________________
446void AliRelAlignerKalmanArray::Browse(TBrowser* b )
447{
448 if (!b) return;
449 if (!fListOfGraphs)
450 {
451 fListOfGraphs=new TList();
452 fListOfGraphs->SetName("graphs");
453 fListOfGraphs->SetOwner(kTRUE);
454 }
455 for (Int_t i=0;i<9;i++)
456 {
457 TGraphErrors* gr = dynamic_cast<TGraphErrors*>(fListOfGraphs->At(i));
458 if (gr) continue;
459 fListOfGraphs->AddAt(MakeGraph(i),i);
460 }
461 b->Add(fListOfGraphs);
462}
463