- Coverity fixes
[u/mrichter/AliRoot.git] / PWGPP / 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
e99fb5c9 209 if(&in == this) return *this;
210 TNamed::operator=(in);
211
51cfc50d 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;
4a84c20d 233}
234
235//______________________________________________________________________________
51cfc50d 236void AliRelAlignerKalmanArray::ClearContents()
4a84c20d 237{
51cfc50d 238 //clear contents of array
239 for (Int_t i=0;i<fSize;i++)
240 {
241 if (fPArray[i]) delete fPArray[i];
242 }
4a84c20d 243}
244
51cfc50d 245//______________________________________________________________________________
246AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
4a84c20d 247{
51cfc50d 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}
4a84c20d 253
51cfc50d 254//______________________________________________________________________________
255AliRelAlignerKalman* 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];
4a84c20d 261}
262
263//______________________________________________________________________________
51cfc50d 264AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
4a84c20d 265{
51cfc50d 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];
4a84c20d 269}
270
271//______________________________________________________________________________
51cfc50d 272AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
4a84c20d 273{
51cfc50d 274 //return aligner in last filled slot
275 if (fSize==0) return NULL;
276 return fPArray[fSize-1];
4a84c20d 277}
278
279//______________________________________________________________________________
51cfc50d 280void AliRelAlignerKalmanArray::Print(Option_t* option) const
4a84c20d 281{
51cfc50d 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 }
4a84c20d 295}
296
297//______________________________________________________________________________
51cfc50d 298Int_t AliRelAlignerKalmanArray::GetEntries() const
4a84c20d 299{
51cfc50d 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;
4a84c20d 308}
309
310//______________________________________________________________________________
51cfc50d 311void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
4a84c20d 312{
51cfc50d 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 }
4a84c20d 321}
322
323//______________________________________________________________________________
51cfc50d 324TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
4a84c20d 325{
51cfc50d 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
63f4d30f 333 Int_t n=GetEntries();
334 TVectorD vx(n);
335 TVectorD vy(n);
336 TVectorD vex(n);
337 TVectorD vey(n);
51cfc50d 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 }
63f4d30f 392 graph->SetName(graphtitle);
51cfc50d 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;
4a84c20d 400}
401
402//______________________________________________________________________________
51cfc50d 403AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
4a84c20d 404{
51cfc50d 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;
4a84c20d 432}
433
434//______________________________________________________________________________
51cfc50d 435void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
4a84c20d 436{
51cfc50d 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;
4a84c20d 446}
447
63f4d30f 448//______________________________________________________________________________
449void 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