1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
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
22 // Origin: Mikolaj Krzewicki, Nikhef, Mikolaj.Krzewicki@cern.ch
24 //////////////////////////////////////////////////////////////////////////////
27 #include <TGraphErrors.h>
31 #include <TObjArray.h>
32 #include <TCollection.h>
33 #include <AliESDEvent.h>
34 #include <AliRelAlignerKalman.h>
35 #include "AliRelAlignerKalmanArray.h"
39 ClassImp(AliRelAlignerKalmanArray)
41 //______________________________________________________________________________
42 AliRelAlignerKalmanArray::AliRelAlignerKalmanArray():
43 TNamed("alignerArray","array of aligners"),
47 fOutRejSigmaOnMerge(10.),
48 fOutRejSigmaOnSmooth(1.),
56 //______________________________________________________________________________
57 AliRelAlignerKalmanArray::AliRelAlignerKalmanArray(Int_t t0, Int_t tend, Int_t slotwidth):
58 TNamed("alignerArray","array of aligners"),
60 fTimebinWidth(slotwidth),
62 fOutRejSigmaOnMerge(10.),
63 fOutRejSigmaOnSmooth(1.),
66 fListOfGraphs(new TList)
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
74 fListOfGraphs->SetName("graphs");
75 fListOfGraphs->SetOwner(kTRUE);
78 //______________________________________________________________________________
79 AliRelAlignerKalmanArray::AliRelAlignerKalmanArray( const AliRelAlignerKalmanArray& in):
80 TNamed(in.GetName(), in.GetTitle()),
82 fTimebinWidth(in.fTimebinWidth),
84 fOutRejSigmaOnMerge(in.fOutRejSigmaOnMerge),
85 fOutRejSigmaOnSmooth(in.fOutRejSigmaOnSmooth),
86 fAlignerTemplate(in.fAlignerTemplate),
88 fListOfGraphs(new TList)
91 fPArray = new AliRelAlignerKalman* [fSize];
92 if (!fPArray) {fSize=0;return;} //if fail
93 for (Int_t i=0;i<fSize;i++)
95 if (in.fPArray[i]) fPArray[i]=new AliRelAlignerKalman(*(in.fPArray[i]));
98 fListOfGraphs->SetName("graphs");
99 fListOfGraphs->SetOwner(kTRUE);
102 //______________________________________________________________________________
103 AliRelAlignerKalmanArray::~AliRelAlignerKalmanArray()
108 delete fListOfGraphs;
111 //______________________________________________________________________________
112 void AliRelAlignerKalmanArray::SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
114 //setup array - clears old content
117 fTimebinWidth=slotwidth;
119 if (slotwidth==0) tmpsize = 1;
120 else tmpsize = (tend-t0)/slotwidth;
124 fPArray=new AliRelAlignerKalman* [tmpsize];
125 if (fPArray) fSize=tmpsize;
128 for (Int_t i=0;i<fSize;i++){fPArray[i]=NULL;}//fill with zeros
131 //______________________________________________________________________________
132 AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAlignerTemplate()
134 //get aligner template
135 return &fAlignerTemplate;
138 //______________________________________________________________________________
139 Long64_t AliRelAlignerKalmanArray::Merge( TCollection* list )
141 //Merge all the arrays
142 //tlihe merge is vertical, meaning matching entries in tree are merged
144 AliRelAlignerKalmanArray *arrayFromList;
147 while ( (arrayFromList = dynamic_cast<AliRelAlignerKalmanArray*>(next())) )
149 if (fT0 != arrayFromList->fT0) continue;
150 if (fTimebinWidth != arrayFromList->fTimebinWidth) continue;
151 if (fSize != arrayFromList->fSize) continue;
153 for (Int_t i=0; i<GetSize(); i++)
155 AliRelAlignerKalman* a1 = fPArray[i];
156 AliRelAlignerKalman* a2 = arrayFromList->fPArray[i];
159 a1->SetRejectOutliers(kFALSE);
160 a1->SetOutRejSigma(fOutRejSigmaOnMerge); a1->Merge(a2);
163 if (!a1 && a2) fPArray[i] = new AliRelAlignerKalman(*a2);
166 fListOfGraphs->Delete();
170 //______________________________________________________________________________
171 Int_t AliRelAlignerKalmanArray::Timebin( UInt_t timestamp ) const
173 //calculate binnumber for given timestamp
174 if (fTimebinWidth==0) return 0;
175 Int_t slot = (timestamp-fT0)/fTimebinWidth; //it's all integers!
179 //______________________________________________________________________________
180 AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(UInt_t ts)
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;
189 fPArray[tb] = new AliRelAlignerKalman( fAlignerTemplate );
190 fPArray[tb]->SetTimeStamp(fT0+tb*fTimebinWidth);
195 //______________________________________________________________________________
196 AliRelAlignerKalman* AliRelAlignerKalmanArray::GetAligner(AliESDEvent* event)
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());
205 //______________________________________________________________________________
206 AliRelAlignerKalmanArray& AliRelAlignerKalmanArray::operator=(const AliRelAlignerKalmanArray& in)
208 //assignment operator
209 if(&in == this) return *this;
210 TNamed::operator=(in);
214 //if sizes different, delete curent and make a new one with new size
217 fPArray = new AliRelAlignerKalman* [in.fSize];
220 fOutRejSigmaOnMerge=in.fOutRejSigmaOnMerge;
221 fOutRejSigmaOnSmooth=in.fOutRejSigmaOnSmooth;
223 if (!fPArray) fSize=0;
224 else fSize = in.fSize;
225 fTimebinWidth = in.fTimebinWidth;
227 for (Int_t i=0;i<fSize;i++)
229 if (in.fPArray[i]) fPArray[i] = new AliRelAlignerKalman(*(in.fPArray[i]));
230 else fPArray[i]=NULL;
235 //______________________________________________________________________________
236 void AliRelAlignerKalmanArray::ClearContents()
238 //clear contents of array
239 for (Int_t i=0;i<fSize;i++)
241 if (fPArray[i]) delete fPArray[i];
245 //______________________________________________________________________________
246 AliRelAlignerKalman* AliRelAlignerKalmanArray::At( Int_t i ) const
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;
254 //______________________________________________________________________________
255 AliRelAlignerKalman* AliRelAlignerKalmanArray::operator[](Int_t i) const
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;
263 //______________________________________________________________________________
264 AliRelAlignerKalman*& AliRelAlignerKalmanArray::operator[](Int_t i)
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];}
271 //______________________________________________________________________________
272 AliRelAlignerKalman* AliRelAlignerKalmanArray::Last() const
274 //return aligner in last filled slot
275 if (fSize==0) return NULL;
276 return fPArray[fSize-1];
279 //______________________________________________________________________________
280 void AliRelAlignerKalmanArray::Print(Option_t* option) const
283 TString optionStr(option);
284 printf( "%s: t0: %i, tend: %i, width: %i, size: %i, entries: %i\n",
286 fT0, (fT0+(fSize)*fTimebinWidth), fTimebinWidth,
287 fSize, GetEntries() );
288 if (optionStr.Contains("a"))
289 for (Int_t i=0; i<fSize; i++)
291 AliRelAlignerKalman* al = fPArray[i];
297 //______________________________________________________________________________
298 Int_t AliRelAlignerKalmanArray::GetEntries() const
300 //get number of filled slots
301 if (!fPArray) return 0;
303 for (Int_t i=0; i<fSize; i++)
305 if (fPArray[i]) entries++;
310 //______________________________________________________________________________
311 void AliRelAlignerKalmanArray::FillTree( TTree* tree ) const
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++)
319 if (al) tree->Fill();
323 //______________________________________________________________________________
324 TGraphErrors* AliRelAlignerKalmanArray::MakeGraph(Int_t iparam) const
326 //return a graph for the iparam-th parameter
327 if (iparam>8 || iparam<0)
329 printf("wrong parameter number. must be from 0-8");
333 Int_t n=GetEntries();
339 for (Int_t i=0; i<fSize; i++)
341 AliRelAlignerKalman* al = fPArray[i];
343 vx(entry) = al->GetTimeStamp();
344 vy(entry) = al->GetStateArr()[iparam];
345 TMatrixDSym* cm = al->GetStateCov();
346 vey(entry) = TMath::Sqrt((*cm)(iparam,iparam));
350 TGraphErrors* graph = new TGraphErrors(vx,vy,vex,vey);
356 graphtitle="rotation \\psi";
357 graphtitley="\\psi [rad]";
360 graphtitle="rotation \\theta";
361 graphtitley="\\theta [rad]";
364 graphtitle="rotation \\phi";
365 graphtitley="\\phi [rad]";
368 graphtitle="shift x";
369 graphtitley="x [cm]";
372 graphtitle="shift y";
373 graphtitley="y [cm]";
376 graphtitle="shift z";
377 graphtitley="z [cm]";
380 graphtitle="TPC vd correction";
381 graphtitley="correction factor []";
384 graphtitle="TPC t0 correction";
385 graphtitley="t0 correction [\\micros]";
388 graphtitle="TPC dv/dy";
389 graphtitley="dv/dy [cm/\\micros/m]";
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);
402 //______________________________________________________________________________
403 AliRelAlignerKalmanArray* AliRelAlignerKalmanArray::MakeSmoothArray() const
405 //return a smoothed version of the data
406 AliRelAlignerKalmanArray* outputarr = new AliRelAlignerKalmanArray(fT0,
407 fT0+fSize*fTimebinWidth,fTimebinWidth);
409 AliRelAlignerKalman* tmpaligner = outputarr->GetAlignerTemplate();
410 tmpaligner->SetOutRejSigma(fOutRejSigmaOnSmooth);
411 //copy the first filled slot
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++)
421 if (!fPArray[i]) continue;
422 PropagateToTime(tmpaligner, fPArray[i]->GetTimeStamp());
423 if (tmpaligner->Merge(fPArray[i]))
424 (*outputarr)[i] = new AliRelAlignerKalman(*tmpaligner);
426 (*outputarr)[i] = new AliRelAlignerKalman(*(fPArray[i]));
429 tmpaligner->Reset(); //clean up
434 //______________________________________________________________________________
435 void AliRelAlignerKalmanArray::PropagateToTime(AliRelAlignerKalman* al, UInt_t timestamp ) const
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
448 //______________________________________________________________________________
449 void AliRelAlignerKalmanArray::Browse(TBrowser* b )
454 fListOfGraphs=new TList();
455 fListOfGraphs->SetName("graphs");
456 fListOfGraphs->SetOwner(kTRUE);
458 for (Int_t i=0;i<9;i++)
460 TGraphErrors* gr = dynamic_cast<TGraphErrors*>(fListOfGraphs->At(i));
462 fListOfGraphs->AddAt(MakeGraph(i),i);
464 b->Add(fListOfGraphs);