Work around when AliFMDCalibZeroSuppression contains
[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"
36
37ClassImp(AliRelAlignerKalmanArray)
38
39//______________________________________________________________________________
40AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
51cfc50d 41 TNamed("alignerArray","array of aligners"),
42 fT0(0),
43 fTimebinWidth(0),
44 fSize(0),
45 fOutRejSigmaOnMerge(10.),
46 fOutRejSigmaOnSmooth(1.),
8478bf3d 47 fAlignerTemplate(),
51cfc50d 48 fPArray(NULL)
4a84c20d 49{
50 //ctor
51}
52
53//______________________________________________________________________________
51cfc50d 54AliRelAlignerKalmanArray::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.),
8478bf3d 61 fAlignerTemplate(),
51cfc50d 62 fPArray(NULL)
4a84c20d 63{
64 //ctor
51cfc50d 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//______________________________________________________________________________
73AliRelAlignerKalmanArray::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),
8478bf3d 80 fAlignerTemplate(in.fAlignerTemplate),
51cfc50d 81 fPArray(NULL)
82{
83 //copy ctor
51cfc50d 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 }
4a84c20d 91}
92
93//______________________________________________________________________________
94AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
95{
96 //dtor
51cfc50d 97 ClearContents();
98 delete [] fPArray;
99}
100
101//______________________________________________________________________________
102void 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//______________________________________________________________________________
122AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAlignerTemplate()
123{
8478bf3d 124 //get aligner template
125 return &fAlignerTemplate;
4a84c20d 126}
127
128//______________________________________________________________________________
129Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
130{
131 //Merge all the arrays
51cfc50d 132 //tlihe merge is vertical, meaning matching entries in tree are merged
4a84c20d 133
134 AliRelAlignerKalmanArray *arrayFromList;
135 if (!list) return 0;
136 TIter next(list);
7371114e 137 while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
4a84c20d 138 {
51cfc50d 139 if (fT0 != arrayFromList->fT0) continue;
140 if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
141 if (fSize != arrayFromList->fSize) continue;
4a84c20d 142
51cfc50d 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 }
4a84c20d 155 }
51cfc50d 156 return 0;
157}
4a84c20d 158
51cfc50d 159//______________________________________________________________________________
160Int_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;
4a84c20d 166}
167
168//______________________________________________________________________________
51cfc50d 169AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
4a84c20d 170{
51cfc50d 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])
4a84c20d 177 {
8478bf3d 178 fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
51cfc50d 179 fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
4a84c20d 180 }
51cfc50d 181 return fPArray[tb];
4a84c20d 182}
183
184//______________________________________________________________________________
51cfc50d 185AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
4a84c20d 186{
51cfc50d 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;
4a84c20d 192}
193
194//______________________________________________________________________________
51cfc50d 195AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
4a84c20d 196{
51cfc50d 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;
4a84c20d 219}
220
221//______________________________________________________________________________
51cfc50d 222void AliRelAlignerKalmanArray::ClearContents()
4a84c20d 223{
51cfc50d 224 //clear contents of array
225 for (Int_t i=0;i<fSize;i++)
226 {
227 if (fPArray[i]) delete fPArray[i];
228 }
4a84c20d 229}
230
51cfc50d 231//______________________________________________________________________________
232AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
4a84c20d 233{
51cfc50d 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}
4a84c20d 239
51cfc50d 240//______________________________________________________________________________
241AliRelAlignerKalman* 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];
4a84c20d 247}
248
249//______________________________________________________________________________
51cfc50d 250AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
4a84c20d 251{
51cfc50d 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];
4a84c20d 255}
256
257//______________________________________________________________________________
51cfc50d 258AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
4a84c20d 259{
51cfc50d 260 //return aligner in last filled slot
261 if (fSize==0) return NULL;
262 return fPArray[fSize-1];
4a84c20d 263}
264
265//______________________________________________________________________________
51cfc50d 266void AliRelAlignerKalmanArray::Print(Option_t* option) const
4a84c20d 267{
51cfc50d 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 }
4a84c20d 281}
282
283//______________________________________________________________________________
51cfc50d 284Int_t AliRelAlignerKalmanArray::GetEntries() const
4a84c20d 285{
51cfc50d 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;
4a84c20d 294}
295
296//______________________________________________________________________________
51cfc50d 297void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
4a84c20d 298{
51cfc50d 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 }
4a84c20d 307}
308
309//______________________________________________________________________________
51cfc50d 310TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
4a84c20d 311{
51cfc50d 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;
4a84c20d 384}
385
386//______________________________________________________________________________
51cfc50d 387AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
4a84c20d 388{
51cfc50d 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;
4a84c20d 416}
417
418//______________________________________________________________________________
51cfc50d 419void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
4a84c20d 420{
51cfc50d 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;
4a84c20d 430}
431