]> 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
1. Increase the angular cut 0.015==> 0.05 to allow bigger fluctuation of drift velocity
[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
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 }
416 char name[1000];
417 sprintf(name,"%s Pad 1D",GetTitle());
418 TH1F * his = new TH1F(name,name,100, min,max);
419 for (Int_t isec = 0; isec < kNsec; isec++) {
4486a91f 420 if (side==1 && isec%36>18) continue;
421 if (side==-1 && isec%36<18) continue;
90127643 422 if (fROC[isec]){
423 for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
424 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
425 for (UInt_t ipad=0; ipad<npads; ipad++){
426 his->Fill(fROC[isec]->GetValue(irow,ipad));
427 }
428 }
429 }
430 }
431 return his;
432}
433
434//_____________________________________________________________________________
200be8a6 435TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
436 //
437 // Make 2D graph
438 // side - specify the side A = 0 C = 1
439 // type - used types of determination of boundaries in z
a6d2bd0c 440 //
200be8a6 441 Float_t kEpsilon = 0.000000000001;
442 TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
443 AliTPCROC * roc = AliTPCROC::Instance();
444 for (Int_t isec=0; isec<72; isec++){
445 if (side==0 && isec%36>=18) continue;
446 if (side>0 && isec%36<18) continue;
447 if (fROC[isec]){
448 AliTPCCalROC * calRoc = fROC[isec];
449 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
450 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
451 if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
452 Float_t xyz[3];
453 roc->GetPositionGlobal(isec,irow,ipad,xyz);
454 Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
455 Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
456 Float_t value = calRoc->GetValue(irow,ipad);
457 his->SetBinContent(binx,biny,value);
458 }
459 }
460 }
461 his->SetXTitle("x (cm)");
462 his->SetYTitle("y (cm)");
463 return his;
464}
465
466
72d0ab7e 467AliTPCCalPad* 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 468 //
469 // Loops over all AliTPCCalROCs and performs a localFit in each ROC
470 // AliTPCCalPad with fit-data is returned
471 // rowRadius and padRadius specifies a window around a given pad in one sector.
472 // The data of this window are fitted with a parabolic function.
473 // This function is evaluated at the pad's position.
474 // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
475 // rowRadius - radius - rows to be used for smoothing
476 // padradius - radius - pads to be used for smoothing
477 // ROCoutlier - map of outliers - pads not to be used for local smoothing
478 // robust - robust method of fitting - (much slower)
479 // chi2Threshold: Threshold for chi2 when EvalRobust is called
480 // robustFraction: Fraction of data that will be used in EvalRobust
481 //
482 //
483 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
484 for (Int_t isec = 0; isec < 72; isec++){
72d0ab7e 485 if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
a6d2bd0c 486 if (PadOutliers)
72d0ab7e 487 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
a6d2bd0c 488 else
72d0ab7e 489 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
a6d2bd0c 490 }
491 return pad;
492}
493
494
6d8681ef 495AliTPCCalPad* 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 496 //
497 // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
498 // AliTPCCalPad with fit-data is returned
499 // chi2Threshold: Threshold for chi2 when EvalRobust is called
500 // robustFraction: Fraction of data that will be used in EvalRobust
501 // chi2Threshold: Threshold for chi2 when EvalRobust is called
502 // robustFraction: Fraction of data that will be used in EvalRobust
6d8681ef 503 // err: error of the data points
504 // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
a6d2bd0c 505 //
506 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
507 TVectorD fitParam(0);
508 TMatrixD covMatrix(0,0);
509 Float_t chi2 = 0;
510 for (Int_t isec = 0; isec < 72; isec++){
511 if (PadOutliers)
6d8681ef 512 GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
a6d2bd0c 513 else
6d8681ef 514 GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
515
516 AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec);
517 pad->SetCalROC(roc);
518 delete roc;
519 if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec);
520 if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec);
a6d2bd0c 521 }
522 return pad;
523}
732e90a8 524//_____________________________________________________________________________
525TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula)
526{
5312f439 527 //
732e90a8 528 // create an array of TFormulas for the each parameter of the fit function
5312f439 529 //
530
531 // split fit string in single parameters
532 // find dimension of the fit:
533 TString fitString(fitFormula);
534 fitString.ReplaceAll("++","#");
535 fitString.ReplaceAll(" ","");
536 TObjArray *arrFitParams = fitString.Tokenize("#");
537 Int_t ndim = arrFitParams->GetEntries();
5312f439 538 //create array of TFormulas to evaluate the parameters
539 TObjArray *arrFitFormulas = new TObjArray(ndim);
540 arrFitFormulas->SetOwner(kTRUE);
541 for (Int_t idim=0;idim<ndim;++idim){
542 TString s=((TObjString*)arrFitParams->At(idim))->GetString();
543 s.ReplaceAll("gx","[0]");
544 s.ReplaceAll("gy","[1]");
545 s.ReplaceAll("lx","[2]");
546 s.ReplaceAll("ly","[3]");
7390f655 547 s.ReplaceAll("sector","[4]");
5312f439 548 arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
549 }
732e90a8 550 delete arrFitParams;
551
552 return arrFitFormulas;
553}
554//_____________________________________________________________________________
555void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results,
556 const Int_t sec, const Int_t row, const Int_t pad)
557{
558 //
559 // evaluate the fit formulas
560 //
561 Int_t ndim=arrFitFormulas.GetEntries();
562 results.ResizeTo(ndim);
563
5312f439 564 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
565 Float_t localXYZ[3];
566 Float_t globalXYZ[3];
732e90a8 567 tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ);
568 tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ);
569 //calculate parameter values
570 for (Int_t idim=0;idim<ndim;++idim){
571 TFormula *f=(TFormula*)arrFitFormulas.At(idim);
572 f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec);
573 results[idim]=f->Eval(0);
574 }
575}
576//_____________________________________________________________________________
577void 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){
578 //
579 // Performs a fit on both sides.
580 // Valid information for the fitFormula are the variables
581 // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName
582 // - sector: the sector number.
583 // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on
584 //
585 // PadOutliers - pads with value !=0 are not used in fitting procedure
586 // chi2Threshold: Threshold for chi2 when EvalRobust is called
587 // robustFraction: Fraction of data that will be used in EvalRobust
588 //
589
590 TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula);
591 Int_t ndim = arrFitFormulas->GetEntries();
592 //resize output data arrays
593 fitParamSideA.ResizeTo(ndim+1);
594 fitParamSideC.ResizeTo(ndim+1);
595 covMatrixSideA.ResizeTo(ndim+1,ndim+1);
596 covMatrixSideC.ResizeTo(ndim+1,ndim+1);
597 // create linear fitter for A- and C- Side
598 TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
599 TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
600 fitterGA->StoreData(kTRUE);
601 fitterGC->StoreData(kTRUE);
602 //parameter values
5312f439 603 TVectorD parValues(ndim);
732e90a8 604
605 AliTPCCalROC *rocErr=0x0;
5312f439 606
607 for (UInt_t isec = 0; isec<kNsec; ++isec){
608 AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
609 AliTPCCalROC *rocData=GetCalROC(isec);
732e90a8 610 if (pointError) rocErr=pointError->GetCalROC(isec);
5312f439 611 if (!rocData) continue;
612 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
613 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
614 //check for outliers
615 if (rocOut && rocOut->GetValue(irow,ipad)) continue;
5312f439 616 //calculate parameter values
732e90a8 617 EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad);
5312f439 618 //get value
619 Float_t value=rocData->GetValue(irow,ipad);
732e90a8 620 //point error
621 Int_t err=1;
622 if (rocErr) {
c018bcd4 623 err=TMath::Nint(rocErr->GetValue(irow,ipad));
732e90a8 624 if (err==0) err=1;
625 }
5312f439 626 //add points to the fitters
627 if (isec/18%2==0){
732e90a8 628 fitterGA->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 629 }else{
732e90a8 630 fitterGC->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 631 }
632 }
633 }
634 }
635 if (robust){
636 fitterGA->EvalRobust(robustFraction);
637 fitterGC->EvalRobust(robustFraction);
638 } else {
639 fitterGA->Eval();
640 fitterGC->Eval();
641 }
642 chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1));
643 chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1));
644 fitterGA->GetParameters(fitParamSideA);
645 fitterGC->GetParameters(fitParamSideC);
646 fitterGA->GetCovarianceMatrix(covMatrixSideA);
647 fitterGC->GetCovarianceMatrix(covMatrixSideC);
648
5312f439 649 delete arrFitFormulas;
650 delete fitterGA;
651 delete fitterGC;
652
653}
732e90a8 654//
655AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC)
656{
657 //
658 //
659 //
660 TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula);
661 Int_t ndim = arrFitFormulas->GetEntries();
662 //check if dimension of fit formula and fit parameters agree
663 if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){
664 printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!");
665 return 0;
666 }
667 //create cal pad
668 AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula));
669 //fill cal pad with fit results if requested
670 for (UInt_t isec = 0; isec<kNsec; ++isec){
671 AliTPCCalROC *roc=pad->GetCalROC(isec);
672 for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) {
673 for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) {
674 const TVectorD *fitPar=0;
675 TVectorD fitResArray;
676 if (isec/18%2==0){
677 fitPar=&fitParamSideA;
678 }else{
679 fitPar=&fitParamSideC;
680 }
681 EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad);
682 for (Int_t idim=0;idim<ndim;++idim)
683 fitResArray(idim)*=(*fitPar)(idim);
684 roc->SetValue(irow,ipad,fitResArray.Sum());
685 }
686 }
687 }
688 delete arrFitFormulas;
689 return pad;
690}
4486a91f 691
692
693
694TCanvas * AliTPCCalPad::MakeReportPadSector(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
a6d2bd0c 695 //
4486a91f 696 // Make a report - cal pads per sector
697 // mean valeus per sector and local X
a6d2bd0c 698 //
4486a91f 699 TH1* his=0;
700 TLegend *legend = 0;
701 TCanvas *canvas = new TCanvas(Form("Sector: %s",varTitle),Form("Sector: %s",varTitle),1500,1100);
702
703 canvas->Divide(2);
704 chain->SetAlias("lX","lx.fElements");
705 //
706 canvas->cd(1);
707 TString strDraw=varName;
708 strDraw+=":lX";
709 legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC A side", varTitle));
710 for (Int_t isec=-1; isec<18; isec+=1){
7d85e147 711 TCut cutSec=Form("sector%%36==%d",isec);
4486a91f 712 cutSec+=cutUser;
713 if (isec==-1) cutSec="sector%36<18";
714 chain->SetMarkerColor(1+(isec+2)%5);
715 chain->SetLineColor(1+(isec+2)%5);
716 chain->SetMarkerStyle(25+(isec+2)%4);
717 //
718 chain->Draw(strDraw.Data(),cutSec,"profgoff");
719 his=(TH1*)chain->GetHistogram()->Clone();
720 delete chain->GetHistogram();
721 his->SetMaximum(max);
722 his->SetMinimum(min);
723 his->GetXaxis()->SetTitle("R (cm)");
724 his->GetYaxis()->SetTitle(axisTitle);
725 his->SetTitle(Form("%s- sector %d",varTitle, isec));
726 his->SetName(Form("%s- sector %d",varTitle, isec));
727 if (isec==-1) his->SetTitle(Form("%s A side",varTitle));
728 if (isec==-1) his->Draw();
729 his->Draw("same");
730 legend->AddEntry(his);
a6d2bd0c 731 }
4486a91f 732 legend->Draw();
733 canvas->cd(2);
734 //
735 legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC C side", varTitle));
736 for (Int_t isec=-1; isec<18; isec+=1){
7d85e147 737 TCut cutSec=Form("(sector+18)%%36==%d",isec);
4486a91f 738 cutSec+=cutUser;
739 if (isec==-1) cutSec="sector%36>18";
740 chain->SetMarkerColor(1+(isec+2)%5);
741 chain->SetLineColor(1+(isec+2)%5);
742 chain->SetMarkerStyle(25+isec%4);
743 //
744 chain->Draw(strDraw.Data(),cutSec,"profgoff");
745 his=(TH1*)chain->GetHistogram()->Clone();
746 delete chain->GetHistogram();
747 his->SetMaximum(max);
748 his->SetMinimum(min);
749 his->GetXaxis()->SetTitle("R (cm)");
750 his->GetYaxis()->SetTitle(axisTitle);
751 his->SetTitle(Form("%s- sector %d",varTitle,isec));
752 his->SetName(Form("%s- sector %d",varTitle,isec));
753 if (isec==-1) his->SetTitle(Form("%s C side",varTitle));
754 if (isec==-1) his->Draw();
755 his->Draw("same");
756 legend->AddEntry(his);
a6d2bd0c 757 }
4486a91f 758 legend->Draw();
759 //
760 //
761 return canvas;
762}
763
764
765TCanvas * AliTPCCalPad::MakeReportPadSector2D(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
766 //
767 // Make a report - cal pads per sector
768 // 2D view
769 // Input tree should be created using AliPreprocesorOnline before
770 //
771 TH1* his=0;
772 TCanvas *canvas = new TCanvas(Form("%s2D",varTitle),Form("%s2D",varTitle),1500,1100);
773 canvas->Divide(2);
774 //
775 TString strDraw=varName;
776 strDraw+=":gy.fElements:gx.fElements>>his(250,-250,250,250,-250,250)";
777 //
778 TVirtualPad * pad=0;
779 pad=canvas->cd(1);
780 pad->SetMargin(0.15,0.15,0.15,0.15);
781 TCut cut=cutUser;
782 chain->Draw(strDraw.Data(),"sector%36<18"+cut,"profgoffcolz2");
783 his=(TH1*)chain->GetHistogram()->Clone();
784 delete chain->GetHistogram();
785 his->SetMaximum(max);
786 his->SetMinimum(min);
787 his->GetXaxis()->SetTitle("x (cm)");
788 his->GetYaxis()->SetTitle("y (cm)");
789 his->GetZaxis()->SetTitle(axisTitle);
790 his->SetTitle(Form("%s A side",varTitle));
791 his->SetName(Form("%s A side",varTitle));
792 his->Draw("colz2");
793 //
794 pad=canvas->cd(2);
795 pad->SetMargin(0.15,0.15,0.15,0.15);
796
797 chain->Draw(strDraw.Data(),"sector%36>=18"+cut,"profgoffcolz2");
798 his=(TH1*)chain->GetHistogram()->Clone();
799 delete chain->GetHistogram();
800 his->SetMaximum(max);
801 his->SetMinimum(min);
802 his->GetXaxis()->SetTitle("x (cm)");
803 his->GetYaxis()->SetTitle("y (cm)");
804 his->GetZaxis()->SetTitle(axisTitle);
805 his->SetTitle(Form("%s C side",varTitle));
806 his->SetName(Form("%s C side",varTitle));
807 his->Draw("colz2");
808 //
809 //
810 return canvas;
a6d2bd0c 811}
586007f3 812
4486a91f 813void AliTPCCalPad::Draw(Option_t* option){
814 //
815 // Draw function - standard 2D view
816 //
817 TH1* his=0;
818 TCanvas *canvas = new TCanvas(Form("%s2D",GetTitle()),Form("%s2D",GetTitle()),900,900);
819 canvas->Divide(2,2);
820 //
821 //
822 TVirtualPad * pad=0;
823 pad=canvas->cd(1);
824 pad->SetMargin(0.15,0.15,0.15,0.15);
825 his=MakeHisto2D(0);
826 his->GetXaxis()->SetTitle("x (cm)");
827 his->GetYaxis()->SetTitle("y (cm)");
828 his->GetZaxis()->SetTitle(GetTitle());
829 his->SetTitle(Form("%s A side",GetTitle()));
830 his->SetName(Form("%s A side",GetTitle()));
831 his->Draw(option);
832 //
833 pad=canvas->cd(2);
834 pad->SetMargin(0.15,0.15,0.15,0.15);
835 his=MakeHisto2D(1);
836 his->GetXaxis()->SetTitle("x (cm)");
837 his->GetYaxis()->SetTitle("y (cm)");
838 his->GetZaxis()->SetTitle(GetTitle());
839 his->SetTitle(Form("%s C side",GetTitle()));
840 his->SetName(Form("%s C side",GetTitle()));
841 his->Draw(option);
842 //
843 pad=canvas->cd(3);
844 pad->SetMargin(0.15,0.15,0.15,0.15);
845 his=MakeHisto1D(-8,8,0,1);
846 his->GetXaxis()->SetTitle(GetTitle());
847 his->SetTitle(Form("%s A side",GetTitle()));
848 his->SetName(Form("%s A side",GetTitle()));
849 his->Draw("err");
850 //
851 pad=canvas->cd(4);
852 pad->SetMargin(0.15,0.15,0.15,0.15);
853 his=MakeHisto1D(-8,8,0,-1);
854 his->GetXaxis()->SetTitle(GetTitle());
855 his->SetTitle(Form("%s C side",GetTitle()));
856 his->SetName(Form("%s C side",GetTitle()));
857 his->Draw("err");
858
859
860}
af6a50bb 861
862
863AliTPCCalPad * AliTPCCalPad::MakeCalPadFromHistoRPHI(TH2 * hisA, TH2* hisC){
864 //
865 // Make cal pad from r-phi histograms
866 //
867 AliTPCROC *proc= AliTPCROC::Instance();
868 AliTPCCalPad *calPad = new AliTPCCalPad("his","his");
869 Float_t globalPos[3];
870 for (Int_t isec=0; isec<72; isec++){
871 AliTPCCalROC* calRoc = calPad->GetCalROC(isec);
872 TH2 * his = ((isec%36<18) ? hisA:hisC);
873 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow+=1){
874 Int_t jrow=irow;
875 if (isec>=36) jrow+=63;
af6a50bb 876 for (UInt_t ipad=0;ipad<proc->GetNPads(isec,irow);ipad+=1){
877 proc->GetPositionGlobal(isec,irow,ipad, globalPos);
878 Double_t phi=TMath::ATan2(globalPos[1],globalPos[0]);
879 //if (phi<0) phi+=TMath::Pi()*2;
880 Int_t bin=his->FindBin(phi,jrow);
881 Float_t value= his->GetBinContent(bin);
882 calRoc->SetValue(irow,ipad,value);
883 }
884 }
885 }
886 return calPad;
887}