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