]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalPad.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Eff C++ warning removal
[u/mrichter/AliRoot.git] / TPC / AliTPCCalPad.cxx
CommitLineData
07627591 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/* $Id$ */
17
18///////////////////////////////////////////////////////////////////////////////
19// //
a6d2bd0c 20// TPC calibration class for parameters which are saved per pad //
21// Each AliTPCCalPad consists of 72 AliTPCCalROC-objects //
07627591 22// //
23///////////////////////////////////////////////////////////////////////////////
24
25#include "AliTPCCalPad.h"
26#include "AliTPCCalROC.h"
184bcc16 27#include <TObjArray.h>
28#include <TAxis.h>
29#include <TGraph.h>
200be8a6 30#include <TGraph2D.h>
31#include <TH2F.h>
586007f3 32#include "TTreeStream.h"
ca5dca67 33#include "TFile.h"
34#include "TKey.h"
5312f439 35#include <TFormula.h>
36#include <TString.h>
37#include <TObjString.h>
72d0ab7e 38#include <iostream>
732e90a8 39#include <AliLog.h>
586007f3 40
07627591 41ClassImp(AliTPCCalPad)
42
43//_____________________________________________________________________________
44AliTPCCalPad::AliTPCCalPad():TNamed()
45{
46 //
47 // AliTPCCalPad default constructor
48 //
49
50 for (Int_t isec = 0; isec < kNsec; isec++) {
51 fROC[isec] = 0;
52 }
53
54}
55
56//_____________________________________________________________________________
57AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title)
58 :TNamed(name,title)
59{
60 //
61 // AliTPCCalPad constructor
62 //
63 for (Int_t isec = 0; isec < kNsec; isec++) {
64 fROC[isec] = new AliTPCCalROC(isec);
65 }
66}
67
68
69//_____________________________________________________________________________
70AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c)
71{
72 //
73 // AliTPCCalPad copy constructor
74 //
75
e6c51e6e 76 for (Int_t isec = 0; isec < kNsec; isec++) {
7fb8d357 77 fROC[isec] = 0;
78 if (c.fROC[isec])
79 fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
e6c51e6e 80 }
07627591 81}
82
184bcc16 83//_____________________________________________________________________________
84AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed()
85{
86 //
87 // AliTPCCalPad default constructor
88 //
89
90 for (Int_t isec = 0; isec < kNsec; isec++) {
91 fROC[isec] = (AliTPCCalROC *)array->At(isec);
92 }
93
94}
95
96
07627591 97///_____________________________________________________________________________
98AliTPCCalPad::~AliTPCCalPad()
99{
100 //
101 // AliTPCCalPad destructor
102 //
103
104 for (Int_t isec = 0; isec < kNsec; isec++) {
105 if (fROC[isec]) {
106 delete fROC[isec];
107 fROC[isec] = 0;
108 }
109 }
110
111}
112
113//_____________________________________________________________________________
114AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c)
115{
116 //
117 // Assignment operator
118 //
119
120 if (this != &c) ((AliTPCCalPad &) c).Copy(*this);
121 return *this;
122
123}
124
125//_____________________________________________________________________________
126void AliTPCCalPad::Copy(TObject &c) const
127{
128 //
129 // Copy function
130 //
131
132 for (Int_t isec = 0; isec < kNsec; isec++) {
133 if (fROC[isec]) {
134 fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]);
135 }
136 }
137 TObject::Copy(c);
138}
184bcc16 139
a6d2bd0c 140
141void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){
142 //
143 // Set AliTPCCalROC copies values from 'roc'
144 // if sector == -1 the sector specified in 'roc' is used
145 // else sector specified in 'roc' is ignored and specified sector is filled
146 //
147 if (sector == -1) sector = roc->GetSector();
72d0ab7e 148 if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector);
a6d2bd0c 149 for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++)
150 fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
151}
152
153
154
90127643 155//_____________________________________________________________________________
156void AliTPCCalPad::Add(Float_t c1)
157{
158 //
a6d2bd0c 159 // add constant c1 to all channels of all ROCs
90127643 160 //
161
162 for (Int_t isec = 0; isec < kNsec; isec++) {
163 if (fROC[isec]){
164 fROC[isec]->Add(c1);
165 }
166 }
167}
168
169//_____________________________________________________________________________
170void AliTPCCalPad::Multiply(Float_t c1)
171{
a6d2bd0c 172 //
173 // multiply each channel of all ROCs with c1
174 //
90127643 175 for (Int_t isec = 0; isec < kNsec; isec++) {
176 if (fROC[isec]){
177 fROC[isec]->Multiply(c1);
178 }
179 }
180}
181
182//_____________________________________________________________________________
183void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
184{
a6d2bd0c 185 //
186 // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
187 // - pad by pad -
188 //
90127643 189 for (Int_t isec = 0; isec < kNsec; isec++) {
79c0f359 190 if (fROC[isec] && pad->GetCalROC(isec)){
90127643 191 fROC[isec]->Add(pad->GetCalROC(isec),c1);
192 }
193 }
194}
195
196//_____________________________________________________________________________
197void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
198{
a6d2bd0c 199 //
200 // multiply each channel of all ROCs with the coresponding channel of 'pad'
201 // - pad by pad -
202 //
203 for (Int_t isec = 0; isec < kNsec; isec++) {
90127643 204 if (fROC[isec]){
205 fROC[isec]->Multiply(pad->GetCalROC(isec));
206 }
207 }
208}
209
210//_____________________________________________________________________________
211void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
212{
a6d2bd0c 213 //
214 // divide each channel of all ROCs by the coresponding channel of 'pad'
215 // - pad by pad -
216 //
90127643 217 for (Int_t isec = 0; isec < kNsec; isec++) {
218 if (fROC[isec]){
219 fROC[isec]->Divide(pad->GetCalROC(isec));
220 }
221 }
222}
223
224//_____________________________________________________________________________
184bcc16 225TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
226 //
227 // type=1 - mean
228 // 2 - median
229 // 3 - LTM
a6d2bd0c 230 //
184bcc16 231 Int_t npoints = 0;
232 for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
233 TGraph * graph = new TGraph(npoints);
234 npoints=0;
235 for (Int_t isec=0;isec<72;isec++){
236 if (!fROC[isec]) continue;
237 if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
238 if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
239 if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
240 npoints++;
241 }
242
243 graph->GetXaxis()->SetTitle("Sector");
244 if (type==0) {
245 graph->GetYaxis()->SetTitle("Mean");
246 graph->SetMarkerStyle(22);
247 }
248 if (type==1) {
249 graph->GetYaxis()->SetTitle("Median");
250 graph->SetMarkerStyle(22);
251 }
252 if (type==2) {
253 graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
254 graph->SetMarkerStyle(24);
255 }
256
257 return graph;
258}
200be8a6 259
90127643 260//_____________________________________________________________________________
261Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
262{
263 //
a6d2bd0c 264 // Calculates mean and RMS of all ROCs
90127643 265 //
266 Double_t sum = 0, sum2 = 0, n=0, val=0;
267 for (Int_t isec = 0; isec < kNsec; isec++) {
268 AliTPCCalROC *calRoc = fROC[isec];
269 if ( calRoc ){
270 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
271 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
272 val = calRoc->GetValue(irow,ipad);
273 sum+=val;
274 sum2+=val*val;
275 n++;
276 }
277 }
278
279 }
280 }
281 Double_t n1 = 1./n;
282 Double_t mean = sum*n1;
283 rms = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
284 return mean;
285}
286
287
288//_____________________________________________________________________________
ca5dca67 289Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
90127643 290{
291 //
292 // return mean of the mean of all ROCs
293 //
294 Double_t arr[kNsec];
295 Int_t n=0;
296 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 297 AliTPCCalROC *calRoc = fROC[isec];
298 if ( calRoc ){
299 AliTPCCalROC* outlierROC = 0;
300 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
301 arr[n] = calRoc->GetMean(outlierROC);
302 n++;
303 }
90127643 304 }
305 return TMath::Mean(n,arr);
306}
307
308//_____________________________________________________________________________
ca5dca67 309Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
90127643 310{
311 //
312 // return mean of the RMS of all ROCs
313 //
314 Double_t arr[kNsec];
315 Int_t n=0;
316 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 317 AliTPCCalROC *calRoc = fROC[isec];
318 if ( calRoc ){
319 AliTPCCalROC* outlierROC = 0;
320 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
321 arr[n] = calRoc->GetRMS(outlierROC);
322 n++;
323 }
90127643 324 }
325 return TMath::Mean(n,arr);
326}
327
328//_____________________________________________________________________________
ca5dca67 329Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
90127643 330{
331 //
332 // return mean of the median of all ROCs
333 //
334 Double_t arr[kNsec];
335 Int_t n=0;
336 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 337 AliTPCCalROC *calRoc = fROC[isec];
338 if ( calRoc ){
339 AliTPCCalROC* outlierROC = 0;
340 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
341 arr[n] = calRoc->GetMedian(outlierROC);
342 n++;
343 }
90127643 344 }
345 return TMath::Mean(n,arr);
346}
347
348//_____________________________________________________________________________
ca5dca67 349Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
90127643 350{
351 //
352 // return mean of the LTM and sigma of all ROCs
353 //
354 Double_t arrm[kNsec];
355 Double_t arrs[kNsec];
356 Double_t *sTemp=0x0;
357 Int_t n=0;
358
359 for (Int_t isec = 0; isec < kNsec; isec++) {
360 AliTPCCalROC *calRoc = fROC[isec];
361 if ( calRoc ){
362 if ( sigma ) sTemp=arrs+n;
ca5dca67 363 AliTPCCalROC* outlierROC = 0;
364 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
365 arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
90127643 366 n++;
367 }
368 }
369 if ( sigma ) *sigma = TMath::Mean(n,arrs);
370 return TMath::Mean(n,arrm);
371}
372
373//_____________________________________________________________________________
374TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){
375 //
376 // make 1D histo
377 // type -1 = user defined range
378 // 0 = nsigma cut nsigma=min
a6d2bd0c 379 //
90127643 380 if (type>=0){
381 if (type==0){
382 // nsigma range
383 Float_t mean = GetMean();
384 Float_t sigma = GetRMS();
385 Float_t nsigma = TMath::Abs(min);
386 min = mean-nsigma*sigma;
387 max = mean+nsigma*sigma;
388 }
389 if (type==1){
390 // fixed range
391 Float_t mean = GetMedian();
392 Float_t delta = min;
393 min = mean-delta;
394 max = mean+delta;
395 }
396 if (type==2){
397 //
398 // LTM mean +- nsigma
399 //
400 Double_t sigma;
401 Float_t mean = GetLTM(&sigma,max);
402 sigma*=min;
403 min = mean-sigma;
404 max = mean+sigma;
405 }
406 }
407 char name[1000];
408 sprintf(name,"%s Pad 1D",GetTitle());
409 TH1F * his = new TH1F(name,name,100, min,max);
410 for (Int_t isec = 0; isec < kNsec; isec++) {
411 if (fROC[isec]){
412 for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
413 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
414 for (UInt_t ipad=0; ipad<npads; ipad++){
415 his->Fill(fROC[isec]->GetValue(irow,ipad));
416 }
417 }
418 }
419 }
420 return his;
421}
422
423//_____________________________________________________________________________
200be8a6 424TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
425 //
426 // Make 2D graph
427 // side - specify the side A = 0 C = 1
428 // type - used types of determination of boundaries in z
a6d2bd0c 429 //
200be8a6 430 Float_t kEpsilon = 0.000000000001;
431 TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
432 AliTPCROC * roc = AliTPCROC::Instance();
433 for (Int_t isec=0; isec<72; isec++){
434 if (side==0 && isec%36>=18) continue;
435 if (side>0 && isec%36<18) continue;
436 if (fROC[isec]){
437 AliTPCCalROC * calRoc = fROC[isec];
438 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
439 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
440 if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
441 Float_t xyz[3];
442 roc->GetPositionGlobal(isec,irow,ipad,xyz);
443 Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
444 Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
445 Float_t value = calRoc->GetValue(irow,ipad);
446 his->SetBinContent(binx,biny,value);
447 }
448 }
449 }
450 his->SetXTitle("x (cm)");
451 his->SetYTitle("y (cm)");
452 return his;
453}
454
455
72d0ab7e 456AliTPCCalPad* AliTPCCalPad::LocalFit(const char* padName, Int_t rowRadius, Int_t padRadius, AliTPCCalPad* PadOutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction, Bool_t printCurrentSector) const {
a6d2bd0c 457 //
458 // Loops over all AliTPCCalROCs and performs a localFit in each ROC
459 // AliTPCCalPad with fit-data is returned
460 // rowRadius and padRadius specifies a window around a given pad in one sector.
461 // The data of this window are fitted with a parabolic function.
462 // This function is evaluated at the pad's position.
463 // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
464 // rowRadius - radius - rows to be used for smoothing
465 // padradius - radius - pads to be used for smoothing
466 // ROCoutlier - map of outliers - pads not to be used for local smoothing
467 // robust - robust method of fitting - (much slower)
468 // chi2Threshold: Threshold for chi2 when EvalRobust is called
469 // robustFraction: Fraction of data that will be used in EvalRobust
470 //
471 //
472 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
473 for (Int_t isec = 0; isec < 72; isec++){
72d0ab7e 474 if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
a6d2bd0c 475 if (PadOutliers)
72d0ab7e 476 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
a6d2bd0c 477 else
72d0ab7e 478 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
a6d2bd0c 479 }
480 return pad;
481}
482
483
6d8681ef 484AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err, TObjArray *fitParArr, TObjArray *fitCovArr){
a6d2bd0c 485 //
486 // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
487 // AliTPCCalPad with fit-data is returned
488 // chi2Threshold: Threshold for chi2 when EvalRobust is called
489 // robustFraction: Fraction of data that will be used in EvalRobust
490 // chi2Threshold: Threshold for chi2 when EvalRobust is called
491 // robustFraction: Fraction of data that will be used in EvalRobust
6d8681ef 492 // err: error of the data points
493 // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
a6d2bd0c 494 //
495 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
496 TVectorD fitParam(0);
497 TMatrixD covMatrix(0,0);
498 Float_t chi2 = 0;
499 for (Int_t isec = 0; isec < 72; isec++){
500 if (PadOutliers)
6d8681ef 501 GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
a6d2bd0c 502 else
6d8681ef 503 GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
504
505 AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec);
506 pad->SetCalROC(roc);
507 delete roc;
508 if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec);
509 if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec);
a6d2bd0c 510 }
511 return pad;
512}
732e90a8 513//_____________________________________________________________________________
514TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula)
515{
5312f439 516 //
732e90a8 517 // create an array of TFormulas for the each parameter of the fit function
5312f439 518 //
519
520 // split fit string in single parameters
521 // find dimension of the fit:
522 TString fitString(fitFormula);
523 fitString.ReplaceAll("++","#");
524 fitString.ReplaceAll(" ","");
525 TObjArray *arrFitParams = fitString.Tokenize("#");
526 Int_t ndim = arrFitParams->GetEntries();
5312f439 527 //create array of TFormulas to evaluate the parameters
528 TObjArray *arrFitFormulas = new TObjArray(ndim);
529 arrFitFormulas->SetOwner(kTRUE);
530 for (Int_t idim=0;idim<ndim;++idim){
531 TString s=((TObjString*)arrFitParams->At(idim))->GetString();
532 s.ReplaceAll("gx","[0]");
533 s.ReplaceAll("gy","[1]");
534 s.ReplaceAll("lx","[2]");
535 s.ReplaceAll("ly","[3]");
7390f655 536 s.ReplaceAll("sector","[4]");
5312f439 537 arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
538 }
732e90a8 539 delete arrFitParams;
540
541 return arrFitFormulas;
542}
543//_____________________________________________________________________________
544void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results,
545 const Int_t sec, const Int_t row, const Int_t pad)
546{
547 //
548 // evaluate the fit formulas
549 //
550 Int_t ndim=arrFitFormulas.GetEntries();
551 results.ResizeTo(ndim);
552
5312f439 553 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
554 Float_t localXYZ[3];
555 Float_t globalXYZ[3];
732e90a8 556 tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ);
557 tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ);
558 //calculate parameter values
559 for (Int_t idim=0;idim<ndim;++idim){
560 TFormula *f=(TFormula*)arrFitFormulas.At(idim);
561 f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec);
562 results[idim]=f->Eval(0);
563 }
564}
565//_____________________________________________________________________________
566void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, const char* fitFormula, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, AliTPCCalPad *pointError, Bool_t robust, Double_t robustFraction){
567 //
568 // Performs a fit on both sides.
569 // Valid information for the fitFormula are the variables
570 // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName
571 // - sector: the sector number.
572 // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on
573 //
574 // PadOutliers - pads with value !=0 are not used in fitting procedure
575 // chi2Threshold: Threshold for chi2 when EvalRobust is called
576 // robustFraction: Fraction of data that will be used in EvalRobust
577 //
578
579 TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula);
580 Int_t ndim = arrFitFormulas->GetEntries();
581 //resize output data arrays
582 fitParamSideA.ResizeTo(ndim+1);
583 fitParamSideC.ResizeTo(ndim+1);
584 covMatrixSideA.ResizeTo(ndim+1,ndim+1);
585 covMatrixSideC.ResizeTo(ndim+1,ndim+1);
586 // create linear fitter for A- and C- Side
587 TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
588 TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
589 fitterGA->StoreData(kTRUE);
590 fitterGC->StoreData(kTRUE);
591 //parameter values
5312f439 592 TVectorD parValues(ndim);
732e90a8 593
594 AliTPCCalROC *rocErr=0x0;
5312f439 595
596 for (UInt_t isec = 0; isec<kNsec; ++isec){
597 AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
598 AliTPCCalROC *rocData=GetCalROC(isec);
732e90a8 599 if (pointError) rocErr=pointError->GetCalROC(isec);
5312f439 600 if (!rocData) continue;
601 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
602 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
603 //check for outliers
604 if (rocOut && rocOut->GetValue(irow,ipad)) continue;
5312f439 605 //calculate parameter values
732e90a8 606 EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad);
5312f439 607 //get value
608 Float_t value=rocData->GetValue(irow,ipad);
732e90a8 609 //point error
610 Int_t err=1;
611 if (rocErr) {
c018bcd4 612 err=TMath::Nint(rocErr->GetValue(irow,ipad));
732e90a8 613 if (err==0) err=1;
614 }
5312f439 615 //add points to the fitters
616 if (isec/18%2==0){
732e90a8 617 fitterGA->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 618 }else{
732e90a8 619 fitterGC->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 620 }
621 }
622 }
623 }
624 if (robust){
625 fitterGA->EvalRobust(robustFraction);
626 fitterGC->EvalRobust(robustFraction);
627 } else {
628 fitterGA->Eval();
629 fitterGC->Eval();
630 }
631 chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1));
632 chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1));
633 fitterGA->GetParameters(fitParamSideA);
634 fitterGC->GetParameters(fitParamSideC);
635 fitterGA->GetCovarianceMatrix(covMatrixSideA);
636 fitterGC->GetCovarianceMatrix(covMatrixSideC);
637
5312f439 638 delete arrFitFormulas;
639 delete fitterGA;
640 delete fitterGC;
641
642}
732e90a8 643//
644AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC)
645{
646 //
647 //
648 //
649 TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula);
650 Int_t ndim = arrFitFormulas->GetEntries();
651 //check if dimension of fit formula and fit parameters agree
652 if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){
653 printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!");
654 return 0;
655 }
656 //create cal pad
657 AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula));
658 //fill cal pad with fit results if requested
659 for (UInt_t isec = 0; isec<kNsec; ++isec){
660 AliTPCCalROC *roc=pad->GetCalROC(isec);
661 for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) {
662 for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) {
663 const TVectorD *fitPar=0;
664 TVectorD fitResArray;
665 if (isec/18%2==0){
666 fitPar=&fitParamSideA;
667 }else{
668 fitPar=&fitParamSideC;
669 }
670 EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad);
671 for (Int_t idim=0;idim<ndim;++idim)
672 fitResArray(idim)*=(*fitPar)(idim);
673 roc->SetValue(irow,ipad,fitResArray.Sum());
674 }
675 }
676 }
677 delete arrFitFormulas;
678 return pad;
679}
5312f439 680/*
a6d2bd0c 681void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, Int_t fitType, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){
682 //
683 // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side
684 // fitType == 0: fit plane function
685 // fitType == 1: fit parabolic function
686 // PadOutliers - pads with value !=0 are not used in fitting procedure
687 // chi2Threshold: Threshold for chi2 when EvalRobust is called
688 // robustFraction: Fraction of data that will be used in EvalRobust
689 //
690 TLinearFitter* fitterGA = 0;
691 TLinearFitter* fitterGC = 0;
692
693 if (fitType == 1) {
694 fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
695 fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
696 }
697 else {
698 fitterGA = new TLinearFitter(3,"x0++x1++x2");
699 fitterGC = new TLinearFitter(3,"x0++x1++x2");
700 }
701 fitterGA->StoreData(kTRUE);
702 fitterGC->StoreData(kTRUE);
703 fitterGA->ClearPoints();
704 fitterGC->ClearPoints();
705 Double_t xx[6];
706 Int_t npointsA=0;
707 Int_t npointsC=0;
708
709 Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position
710 Float_t lx, ly; // pads position
711
712 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
713
714 // loop over all sectors and pads and read data into fitterGA and fitterGC
715 if (fitType == 1) {
716 // parabolic fit
717 fitParamSideA.ResizeTo(6);
718 fitParamSideC.ResizeTo(6);
719 covMatrixSideA.ResizeTo(6,6);
720 covMatrixSideC.ResizeTo(6,6);
721 for (UInt_t isec = 0; isec<72; isec++){
722 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
723 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
724 // fill fitterG
725 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
726 lx = localXY[0];
727 ly = localXY[1];
728 xx[0] = 1;
729 xx[1] = lx;
730 xx[2] = ly;
731 xx[3] = lx*lx;
732 xx[4] = ly*ly;
733 xx[5] = lx*ly;
734 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
735 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
736// sector 0 - 17: IROC, A
737// sector 18 - 35: IROC, C
738// sector 36 - 53: OROC, A
739// sector 54 - 71: CROC, C
740 if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A
741 npointsA++;
742 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
743 }
744 else { // side C
745 npointsC++;
746 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
747 }
748 }
749 }
750 }
751 }
752 }
753 else {
754 // linear fit
755 fitParamSideA.ResizeTo(3);
756 fitParamSideC.ResizeTo(3);
757 covMatrixSideA.ResizeTo(3,3);
758 covMatrixSideC.ResizeTo(3,3);
759
760 for (UInt_t isec = 0; isec<72; isec++){
761 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
762 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
763 // fill fitterG
764 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
765 lx = localXY[0];
766 ly = localXY[1];
767 xx[0] = 1;
768 xx[1] = lx;
769 xx[2] = ly;
770 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
771 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
772// sector 0 - 17: IROC, A
773// sector 18 - 35: IROC, C
774// sector 36 - 53: OROC, A
775// sector 54 - 71: CROC, C
776 if (isec <= 17 || (isec >= 36 && isec <= 53)) {
777 // Side A
778 npointsA++;
779 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
780 }
781 else {
782 // side C
783 npointsC++;
784 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
785 }
786 }
787 }
788 }
789 }
790 }
791
792 fitterGA->Eval();
793 fitterGC->Eval();
794 fitterGA->GetParameters(fitParamSideA);
795 fitterGC->GetParameters(fitParamSideC);
796 fitterGA->GetCovarianceMatrix(covMatrixSideA);
797 fitterGC->GetCovarianceMatrix(covMatrixSideC);
798 if (fitType == 1){
799 chi2SideA = fitterGA->GetChisquare()/(npointsA-6.);
800 chi2SideC = fitterGC->GetChisquare()/(npointsC-6.);
801 }
802 else {
803 chi2SideA = fitterGA->GetChisquare()/(npointsA-3.);
804 chi2SideC = fitterGC->GetChisquare()/(npointsC-3.);
805 }
806 if (robust && chi2SideA > chi2Threshold) {
807 // std::cout << "robust fitter called... " << std::endl;
808 fitterGA->EvalRobust(robustFraction);
809 fitterGA->GetParameters(fitParamSideA);
810 }
811 if (robust && chi2SideC > chi2Threshold) {
812 // std::cout << "robust fitter called... " << std::endl;
813 fitterGC->EvalRobust(robustFraction);
814 fitterGC->GetParameters(fitParamSideC);
815 }
816 delete fitterGA;
817 delete fitterGC;
818}
5312f439 819*/
586007f3 820