]> 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 add Reset function to CalPad and CalROC o Add functionality to AliTPCdataQA - Reset...
[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>
0e00d630 48#include "AliTPCPreprocessorOnline.h"
49#include "AliTPCCalibViewer.h"
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//_____________________________________________________________________________
f345579a 93AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed(array->GetName(),array->GetName())
184bcc16 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
99d61d46 162Bool_t AliTPCCalPad::MedianFilter(Int_t deltaRow, Int_t deltaPad, AliTPCCalPad*outlierPad, Bool_t doEdge){
0e00d630 163 //
164 // replace constent with median in the neigborhood
165 //
166 Bool_t isOK=kTRUE;
167 for (Int_t isec = 0; isec < kNsec; isec++) {
99d61d46 168 AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
0e00d630 169 if (fROC[isec]){
99d61d46 170 isOK&=fROC[isec]->MedianFilter(deltaRow,deltaPad,outlierROC,doEdge);
0e00d630 171 }
172 }
173 return isOK;
174}
a6d2bd0c 175
99d61d46 176Bool_t AliTPCCalPad::LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type, AliTPCCalPad*outlierPad, Bool_t doEdge){
0e00d630 177 //
178 // replace constent with LTM statistic in neigborhood
179 //
180 Bool_t isOK=kTRUE;
181 for (Int_t isec = 0; isec < kNsec; isec++) {
99d61d46 182 AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
0e00d630 183 if (fROC[isec]){
99d61d46 184 isOK&=fROC[isec]->LTMFilter(deltaRow, deltaPad,fraction,type,outlierROC,doEdge);
0e00d630 185 }
186 }
187 return isOK;
188}
a6d2bd0c 189
f787f642 190Bool_t AliTPCCalPad::Convolute(Double_t sigmaPad, Double_t sigmaRow, AliTPCCalPad*outlierPad, TF1 *fpad, TF1 *frow){
191 //
192 // replace constent with median in the neigborhood
193 //
194 Bool_t isOK=kTRUE;
195 for (Int_t isec = 0; isec < kNsec; isec++) {
196 AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
197 if (fROC[isec]){
198 isOK&=fROC[isec]->Convolute(sigmaPad,sigmaRow,outlierROC,fpad,frow);
199 }
200 }
201 return isOK;
202}
203
204
90127643 205//_____________________________________________________________________________
206void AliTPCCalPad::Add(Float_t c1)
207{
208 //
a6d2bd0c 209 // add constant c1 to all channels of all ROCs
90127643 210 //
211
212 for (Int_t isec = 0; isec < kNsec; isec++) {
213 if (fROC[isec]){
214 fROC[isec]->Add(c1);
215 }
216 }
217}
218
219//_____________________________________________________________________________
220void AliTPCCalPad::Multiply(Float_t c1)
221{
a6d2bd0c 222 //
223 // multiply each channel of all ROCs with c1
224 //
90127643 225 for (Int_t isec = 0; isec < kNsec; isec++) {
226 if (fROC[isec]){
227 fROC[isec]->Multiply(c1);
228 }
229 }
230}
231
232//_____________________________________________________________________________
233void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
234{
a6d2bd0c 235 //
236 // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
237 // - pad by pad -
238 //
90127643 239 for (Int_t isec = 0; isec < kNsec; isec++) {
79c0f359 240 if (fROC[isec] && pad->GetCalROC(isec)){
90127643 241 fROC[isec]->Add(pad->GetCalROC(isec),c1);
242 }
243 }
244}
245
246//_____________________________________________________________________________
247void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
248{
a6d2bd0c 249 //
250 // multiply each channel of all ROCs with the coresponding channel of 'pad'
251 // - pad by pad -
252 //
253 for (Int_t isec = 0; isec < kNsec; isec++) {
90127643 254 if (fROC[isec]){
255 fROC[isec]->Multiply(pad->GetCalROC(isec));
256 }
257 }
258}
259
260//_____________________________________________________________________________
261void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
262{
a6d2bd0c 263 //
264 // divide each channel of all ROCs by the coresponding channel of 'pad'
265 // - pad by pad -
266 //
90127643 267 for (Int_t isec = 0; isec < kNsec; isec++) {
268 if (fROC[isec]){
269 fROC[isec]->Divide(pad->GetCalROC(isec));
270 }
271 }
272}
273
f345579a 274//_____________________________________________________________________________
275void AliTPCCalPad::Reset()
276{
277 //
278 // Reset all cal Rocs
279 //
280 for (Int_t isec = 0; isec < kNsec; isec++) {
281 if (fROC[isec]){
282 fROC[isec]->Reset();
283 }
284 }
285}
286
90127643 287//_____________________________________________________________________________
184bcc16 288TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
289 //
290 // type=1 - mean
291 // 2 - median
292 // 3 - LTM
a6d2bd0c 293 //
184bcc16 294 Int_t npoints = 0;
295 for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
296 TGraph * graph = new TGraph(npoints);
297 npoints=0;
298 for (Int_t isec=0;isec<72;isec++){
299 if (!fROC[isec]) continue;
300 if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
301 if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
302 if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
303 npoints++;
304 }
305
306 graph->GetXaxis()->SetTitle("Sector");
307 if (type==0) {
308 graph->GetYaxis()->SetTitle("Mean");
309 graph->SetMarkerStyle(22);
310 }
311 if (type==1) {
312 graph->GetYaxis()->SetTitle("Median");
313 graph->SetMarkerStyle(22);
314 }
315 if (type==2) {
316 graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
317 graph->SetMarkerStyle(24);
318 }
319
320 return graph;
321}
200be8a6 322
90127643 323//_____________________________________________________________________________
324Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
325{
326 //
a6d2bd0c 327 // Calculates mean and RMS of all ROCs
90127643 328 //
329 Double_t sum = 0, sum2 = 0, n=0, val=0;
330 for (Int_t isec = 0; isec < kNsec; isec++) {
331 AliTPCCalROC *calRoc = fROC[isec];
332 if ( calRoc ){
333 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
334 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
335 val = calRoc->GetValue(irow,ipad);
336 sum+=val;
337 sum2+=val*val;
338 n++;
339 }
340 }
341
342 }
343 }
344 Double_t n1 = 1./n;
345 Double_t mean = sum*n1;
346 rms = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
347 return mean;
348}
349
350
351//_____________________________________________________________________________
ca5dca67 352Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
90127643 353{
354 //
355 // return mean of the mean of all ROCs
356 //
357 Double_t arr[kNsec];
358 Int_t n=0;
359 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 360 AliTPCCalROC *calRoc = fROC[isec];
361 if ( calRoc ){
362 AliTPCCalROC* outlierROC = 0;
363 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
364 arr[n] = calRoc->GetMean(outlierROC);
365 n++;
366 }
90127643 367 }
368 return TMath::Mean(n,arr);
369}
370
371//_____________________________________________________________________________
ca5dca67 372Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
90127643 373{
374 //
375 // return mean of the RMS of all ROCs
376 //
377 Double_t arr[kNsec];
378 Int_t n=0;
379 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 380 AliTPCCalROC *calRoc = fROC[isec];
381 if ( calRoc ){
382 AliTPCCalROC* outlierROC = 0;
383 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
384 arr[n] = calRoc->GetRMS(outlierROC);
385 n++;
386 }
90127643 387 }
388 return TMath::Mean(n,arr);
389}
390
391//_____________________________________________________________________________
ca5dca67 392Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
90127643 393{
394 //
395 // return mean of the median of all ROCs
396 //
397 Double_t arr[kNsec];
398 Int_t n=0;
399 for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67 400 AliTPCCalROC *calRoc = fROC[isec];
401 if ( calRoc ){
402 AliTPCCalROC* outlierROC = 0;
403 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
404 arr[n] = calRoc->GetMedian(outlierROC);
405 n++;
406 }
90127643 407 }
408 return TMath::Mean(n,arr);
409}
410
411//_____________________________________________________________________________
ca5dca67 412Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
90127643 413{
414 //
415 // return mean of the LTM and sigma of all ROCs
416 //
417 Double_t arrm[kNsec];
418 Double_t arrs[kNsec];
419 Double_t *sTemp=0x0;
420 Int_t n=0;
421
422 for (Int_t isec = 0; isec < kNsec; isec++) {
423 AliTPCCalROC *calRoc = fROC[isec];
424 if ( calRoc ){
425 if ( sigma ) sTemp=arrs+n;
ca5dca67 426 AliTPCCalROC* outlierROC = 0;
427 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
428 arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
90127643 429 n++;
430 }
431 }
432 if ( sigma ) *sigma = TMath::Mean(n,arrs);
433 return TMath::Mean(n,arrm);
434}
435
436//_____________________________________________________________________________
4486a91f 437TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type, Int_t side){
90127643 438 //
439 // make 1D histo
440 // type -1 = user defined range
441 // 0 = nsigma cut nsigma=min
a6d2bd0c 442 //
90127643 443 if (type>=0){
444 if (type==0){
445 // nsigma range
446 Float_t mean = GetMean();
447 Float_t sigma = GetRMS();
448 Float_t nsigma = TMath::Abs(min);
449 min = mean-nsigma*sigma;
450 max = mean+nsigma*sigma;
451 }
452 if (type==1){
453 // fixed range
454 Float_t mean = GetMedian();
455 Float_t delta = min;
456 min = mean-delta;
457 max = mean+delta;
458 }
459 if (type==2){
460 //
461 // LTM mean +- nsigma
462 //
463 Double_t sigma;
464 Float_t mean = GetLTM(&sigma,max);
465 sigma*=min;
466 min = mean-sigma;
467 max = mean+sigma;
468 }
469 }
5dbad769 470 TString name=Form("%s Pad 1D",GetTitle());
471 TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
90127643 472 for (Int_t isec = 0; isec < kNsec; isec++) {
4486a91f 473 if (side==1 && isec%36>18) continue;
474 if (side==-1 && isec%36<18) continue;
90127643 475 if (fROC[isec]){
476 for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
477 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
478 for (UInt_t ipad=0; ipad<npads; ipad++){
479 his->Fill(fROC[isec]->GetValue(irow,ipad));
480 }
481 }
482 }
483 }
484 return his;
485}
486
487//_____________________________________________________________________________
200be8a6 488TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
489 //
490 // Make 2D graph
491 // side - specify the side A = 0 C = 1
492 // type - used types of determination of boundaries in z
a6d2bd0c 493 //
200be8a6 494 Float_t kEpsilon = 0.000000000001;
495 TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
496 AliTPCROC * roc = AliTPCROC::Instance();
497 for (Int_t isec=0; isec<72; isec++){
498 if (side==0 && isec%36>=18) continue;
499 if (side>0 && isec%36<18) continue;
500 if (fROC[isec]){
501 AliTPCCalROC * calRoc = fROC[isec];
502 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
503 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
504 if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
505 Float_t xyz[3];
506 roc->GetPositionGlobal(isec,irow,ipad,xyz);
507 Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
508 Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
509 Float_t value = calRoc->GetValue(irow,ipad);
510 his->SetBinContent(binx,biny,value);
511 }
512 }
513 }
514 his->SetXTitle("x (cm)");
515 his->SetYTitle("y (cm)");
516 return his;
517}
518
519
72d0ab7e 520AliTPCCalPad* 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 521 //
522 // Loops over all AliTPCCalROCs and performs a localFit in each ROC
523 // AliTPCCalPad with fit-data is returned
524 // rowRadius and padRadius specifies a window around a given pad in one sector.
525 // The data of this window are fitted with a parabolic function.
526 // This function is evaluated at the pad's position.
527 // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
528 // rowRadius - radius - rows to be used for smoothing
529 // padradius - radius - pads to be used for smoothing
530 // ROCoutlier - map of outliers - pads not to be used for local smoothing
531 // robust - robust method of fitting - (much slower)
532 // chi2Threshold: Threshold for chi2 when EvalRobust is called
533 // robustFraction: Fraction of data that will be used in EvalRobust
534 //
535 //
536 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
537 for (Int_t isec = 0; isec < 72; isec++){
72d0ab7e 538 if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
a6d2bd0c 539 if (PadOutliers)
72d0ab7e 540 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
a6d2bd0c 541 else
72d0ab7e 542 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
a6d2bd0c 543 }
544 return pad;
545}
546
547
6d8681ef 548AliTPCCalPad* 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 549 //
550 // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
551 // AliTPCCalPad with fit-data is returned
552 // chi2Threshold: Threshold for chi2 when EvalRobust is called
553 // robustFraction: Fraction of data that will be used in EvalRobust
554 // chi2Threshold: Threshold for chi2 when EvalRobust is called
555 // robustFraction: Fraction of data that will be used in EvalRobust
6d8681ef 556 // err: error of the data points
557 // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
a6d2bd0c 558 //
559 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
560 TVectorD fitParam(0);
561 TMatrixD covMatrix(0,0);
562 Float_t chi2 = 0;
563 for (Int_t isec = 0; isec < 72; isec++){
564 if (PadOutliers)
6d8681ef 565 GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
a6d2bd0c 566 else
6d8681ef 567 GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
568
569 AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec);
570 pad->SetCalROC(roc);
571 delete roc;
572 if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec);
573 if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec);
a6d2bd0c 574 }
575 return pad;
576}
732e90a8 577//_____________________________________________________________________________
578TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula)
579{
5312f439 580 //
732e90a8 581 // create an array of TFormulas for the each parameter of the fit function
5312f439 582 //
583
584 // split fit string in single parameters
585 // find dimension of the fit:
586 TString fitString(fitFormula);
587 fitString.ReplaceAll("++","#");
588 fitString.ReplaceAll(" ","");
589 TObjArray *arrFitParams = fitString.Tokenize("#");
590 Int_t ndim = arrFitParams->GetEntries();
5312f439 591 //create array of TFormulas to evaluate the parameters
592 TObjArray *arrFitFormulas = new TObjArray(ndim);
593 arrFitFormulas->SetOwner(kTRUE);
594 for (Int_t idim=0;idim<ndim;++idim){
595 TString s=((TObjString*)arrFitParams->At(idim))->GetString();
596 s.ReplaceAll("gx","[0]");
597 s.ReplaceAll("gy","[1]");
598 s.ReplaceAll("lx","[2]");
599 s.ReplaceAll("ly","[3]");
7390f655 600 s.ReplaceAll("sector","[4]");
5312f439 601 arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
602 }
732e90a8 603 delete arrFitParams;
604
605 return arrFitFormulas;
606}
607//_____________________________________________________________________________
608void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results,
609 const Int_t sec, const Int_t row, const Int_t pad)
610{
611 //
612 // evaluate the fit formulas
613 //
614 Int_t ndim=arrFitFormulas.GetEntries();
615 results.ResizeTo(ndim);
616
5312f439 617 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
618 Float_t localXYZ[3];
619 Float_t globalXYZ[3];
732e90a8 620 tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ);
621 tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ);
622 //calculate parameter values
623 for (Int_t idim=0;idim<ndim;++idim){
624 TFormula *f=(TFormula*)arrFitFormulas.At(idim);
625 f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec);
626 results[idim]=f->Eval(0);
627 }
628}
629//_____________________________________________________________________________
630void 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){
631 //
632 // Performs a fit on both sides.
633 // Valid information for the fitFormula are the variables
634 // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName
635 // - sector: the sector number.
636 // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on
637 //
638 // PadOutliers - pads with value !=0 are not used in fitting procedure
639 // chi2Threshold: Threshold for chi2 when EvalRobust is called
640 // robustFraction: Fraction of data that will be used in EvalRobust
641 //
642
643 TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula);
644 Int_t ndim = arrFitFormulas->GetEntries();
645 //resize output data arrays
646 fitParamSideA.ResizeTo(ndim+1);
647 fitParamSideC.ResizeTo(ndim+1);
648 covMatrixSideA.ResizeTo(ndim+1,ndim+1);
649 covMatrixSideC.ResizeTo(ndim+1,ndim+1);
650 // create linear fitter for A- and C- Side
651 TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
652 TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
653 fitterGA->StoreData(kTRUE);
654 fitterGC->StoreData(kTRUE);
655 //parameter values
5312f439 656 TVectorD parValues(ndim);
732e90a8 657
658 AliTPCCalROC *rocErr=0x0;
5312f439 659
660 for (UInt_t isec = 0; isec<kNsec; ++isec){
661 AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
662 AliTPCCalROC *rocData=GetCalROC(isec);
732e90a8 663 if (pointError) rocErr=pointError->GetCalROC(isec);
5312f439 664 if (!rocData) continue;
665 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
666 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
667 //check for outliers
668 if (rocOut && rocOut->GetValue(irow,ipad)) continue;
5312f439 669 //calculate parameter values
732e90a8 670 EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad);
5312f439 671 //get value
672 Float_t value=rocData->GetValue(irow,ipad);
732e90a8 673 //point error
674 Int_t err=1;
675 if (rocErr) {
c018bcd4 676 err=TMath::Nint(rocErr->GetValue(irow,ipad));
732e90a8 677 if (err==0) err=1;
678 }
5312f439 679 //add points to the fitters
680 if (isec/18%2==0){
732e90a8 681 fitterGA->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 682 }else{
732e90a8 683 fitterGC->AddPoint(parValues.GetMatrixArray(),value,err);
5312f439 684 }
685 }
686 }
687 }
688 if (robust){
689 fitterGA->EvalRobust(robustFraction);
690 fitterGC->EvalRobust(robustFraction);
691 } else {
692 fitterGA->Eval();
693 fitterGC->Eval();
694 }
695 chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1));
696 chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1));
697 fitterGA->GetParameters(fitParamSideA);
698 fitterGC->GetParameters(fitParamSideC);
699 fitterGA->GetCovarianceMatrix(covMatrixSideA);
700 fitterGC->GetCovarianceMatrix(covMatrixSideC);
701
5312f439 702 delete arrFitFormulas;
703 delete fitterGA;
704 delete fitterGC;
705
706}
732e90a8 707//
708AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC)
709{
710 //
711 //
712 //
713 TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula);
714 Int_t ndim = arrFitFormulas->GetEntries();
715 //check if dimension of fit formula and fit parameters agree
716 if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){
717 printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!");
718 return 0;
719 }
720 //create cal pad
721 AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula));
722 //fill cal pad with fit results if requested
723 for (UInt_t isec = 0; isec<kNsec; ++isec){
724 AliTPCCalROC *roc=pad->GetCalROC(isec);
725 for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) {
726 for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) {
727 const TVectorD *fitPar=0;
728 TVectorD fitResArray;
729 if (isec/18%2==0){
730 fitPar=&fitParamSideA;
731 }else{
732 fitPar=&fitParamSideC;
733 }
734 EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad);
735 for (Int_t idim=0;idim<ndim;++idim)
736 fitResArray(idim)*=(*fitPar)(idim);
737 roc->SetValue(irow,ipad,fitResArray.Sum());
738 }
739 }
740 }
741 delete arrFitFormulas;
742 return pad;
743}
4486a91f 744
745
746
747TCanvas * AliTPCCalPad::MakeReportPadSector(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
a6d2bd0c 748 //
4486a91f 749 // Make a report - cal pads per sector
750 // mean valeus per sector and local X
a6d2bd0c 751 //
4486a91f 752 TH1* his=0;
753 TLegend *legend = 0;
754 TCanvas *canvas = new TCanvas(Form("Sector: %s",varTitle),Form("Sector: %s",varTitle),1500,1100);
755
756 canvas->Divide(2);
757 chain->SetAlias("lX","lx.fElements");
758 //
759 canvas->cd(1);
760 TString strDraw=varName;
761 strDraw+=":lX";
762 legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC A side", varTitle));
763 for (Int_t isec=-1; isec<18; isec+=1){
7d85e147 764 TCut cutSec=Form("sector%%36==%d",isec);
4486a91f 765 cutSec+=cutUser;
766 if (isec==-1) cutSec="sector%36<18";
767 chain->SetMarkerColor(1+(isec+2)%5);
768 chain->SetLineColor(1+(isec+2)%5);
769 chain->SetMarkerStyle(25+(isec+2)%4);
770 //
771 chain->Draw(strDraw.Data(),cutSec,"profgoff");
772 his=(TH1*)chain->GetHistogram()->Clone();
773 delete chain->GetHistogram();
774 his->SetMaximum(max);
775 his->SetMinimum(min);
776 his->GetXaxis()->SetTitle("R (cm)");
777 his->GetYaxis()->SetTitle(axisTitle);
778 his->SetTitle(Form("%s- sector %d",varTitle, isec));
779 his->SetName(Form("%s- sector %d",varTitle, isec));
780 if (isec==-1) his->SetTitle(Form("%s A side",varTitle));
781 if (isec==-1) his->Draw();
782 his->Draw("same");
783 legend->AddEntry(his);
a6d2bd0c 784 }
4486a91f 785 legend->Draw();
786 canvas->cd(2);
787 //
788 legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC C side", varTitle));
789 for (Int_t isec=-1; isec<18; isec+=1){
7d85e147 790 TCut cutSec=Form("(sector+18)%%36==%d",isec);
4486a91f 791 cutSec+=cutUser;
792 if (isec==-1) cutSec="sector%36>18";
793 chain->SetMarkerColor(1+(isec+2)%5);
794 chain->SetLineColor(1+(isec+2)%5);
795 chain->SetMarkerStyle(25+isec%4);
796 //
797 chain->Draw(strDraw.Data(),cutSec,"profgoff");
798 his=(TH1*)chain->GetHistogram()->Clone();
799 delete chain->GetHistogram();
800 his->SetMaximum(max);
801 his->SetMinimum(min);
802 his->GetXaxis()->SetTitle("R (cm)");
803 his->GetYaxis()->SetTitle(axisTitle);
804 his->SetTitle(Form("%s- sector %d",varTitle,isec));
805 his->SetName(Form("%s- sector %d",varTitle,isec));
806 if (isec==-1) his->SetTitle(Form("%s C side",varTitle));
807 if (isec==-1) his->Draw();
808 his->Draw("same");
809 legend->AddEntry(his);
a6d2bd0c 810 }
4486a91f 811 legend->Draw();
812 //
813 //
814 return canvas;
815}
816
817
818TCanvas * AliTPCCalPad::MakeReportPadSector2D(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
819 //
820 // Make a report - cal pads per sector
821 // 2D view
822 // Input tree should be created using AliPreprocesorOnline before
823 //
824 TH1* his=0;
825 TCanvas *canvas = new TCanvas(Form("%s2D",varTitle),Form("%s2D",varTitle),1500,1100);
826 canvas->Divide(2);
827 //
828 TString strDraw=varName;
829 strDraw+=":gy.fElements:gx.fElements>>his(250,-250,250,250,-250,250)";
830 //
831 TVirtualPad * pad=0;
832 pad=canvas->cd(1);
833 pad->SetMargin(0.15,0.15,0.15,0.15);
834 TCut cut=cutUser;
835 chain->Draw(strDraw.Data(),"sector%36<18"+cut,"profgoffcolz2");
836 his=(TH1*)chain->GetHistogram()->Clone();
837 delete chain->GetHistogram();
838 his->SetMaximum(max);
839 his->SetMinimum(min);
840 his->GetXaxis()->SetTitle("x (cm)");
841 his->GetYaxis()->SetTitle("y (cm)");
842 his->GetZaxis()->SetTitle(axisTitle);
843 his->SetTitle(Form("%s A side",varTitle));
844 his->SetName(Form("%s A side",varTitle));
845 his->Draw("colz2");
846 //
847 pad=canvas->cd(2);
848 pad->SetMargin(0.15,0.15,0.15,0.15);
849
850 chain->Draw(strDraw.Data(),"sector%36>=18"+cut,"profgoffcolz2");
851 his=(TH1*)chain->GetHistogram()->Clone();
852 delete chain->GetHistogram();
853 his->SetMaximum(max);
854 his->SetMinimum(min);
855 his->GetXaxis()->SetTitle("x (cm)");
856 his->GetYaxis()->SetTitle("y (cm)");
857 his->GetZaxis()->SetTitle(axisTitle);
858 his->SetTitle(Form("%s C side",varTitle));
859 his->SetName(Form("%s C side",varTitle));
860 his->Draw("colz2");
861 //
862 //
863 return canvas;
a6d2bd0c 864}
586007f3 865
4486a91f 866void AliTPCCalPad::Draw(Option_t* option){
867 //
868 // Draw function - standard 2D view
869 //
870 TH1* his=0;
871 TCanvas *canvas = new TCanvas(Form("%s2D",GetTitle()),Form("%s2D",GetTitle()),900,900);
872 canvas->Divide(2,2);
873 //
874 //
875 TVirtualPad * pad=0;
876 pad=canvas->cd(1);
877 pad->SetMargin(0.15,0.15,0.15,0.15);
878 his=MakeHisto2D(0);
879 his->GetXaxis()->SetTitle("x (cm)");
880 his->GetYaxis()->SetTitle("y (cm)");
881 his->GetZaxis()->SetTitle(GetTitle());
882 his->SetTitle(Form("%s A side",GetTitle()));
883 his->SetName(Form("%s A side",GetTitle()));
884 his->Draw(option);
885 //
886 pad=canvas->cd(2);
887 pad->SetMargin(0.15,0.15,0.15,0.15);
888 his=MakeHisto2D(1);
889 his->GetXaxis()->SetTitle("x (cm)");
890 his->GetYaxis()->SetTitle("y (cm)");
891 his->GetZaxis()->SetTitle(GetTitle());
892 his->SetTitle(Form("%s C side",GetTitle()));
893 his->SetName(Form("%s C side",GetTitle()));
894 his->Draw(option);
895 //
896 pad=canvas->cd(3);
897 pad->SetMargin(0.15,0.15,0.15,0.15);
898 his=MakeHisto1D(-8,8,0,1);
899 his->GetXaxis()->SetTitle(GetTitle());
900 his->SetTitle(Form("%s A side",GetTitle()));
901 his->SetName(Form("%s A side",GetTitle()));
902 his->Draw("err");
903 //
904 pad=canvas->cd(4);
905 pad->SetMargin(0.15,0.15,0.15,0.15);
906 his=MakeHisto1D(-8,8,0,-1);
907 his->GetXaxis()->SetTitle(GetTitle());
908 his->SetTitle(Form("%s C side",GetTitle()));
909 his->SetName(Form("%s C side",GetTitle()));
910 his->Draw("err");
911
912
913}
af6a50bb 914
915
916AliTPCCalPad * AliTPCCalPad::MakeCalPadFromHistoRPHI(TH2 * hisA, TH2* hisC){
917 //
918 // Make cal pad from r-phi histograms
919 //
920 AliTPCROC *proc= AliTPCROC::Instance();
921 AliTPCCalPad *calPad = new AliTPCCalPad("his","his");
922 Float_t globalPos[3];
923 for (Int_t isec=0; isec<72; isec++){
924 AliTPCCalROC* calRoc = calPad->GetCalROC(isec);
925 TH2 * his = ((isec%36<18) ? hisA:hisC);
926 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow+=1){
927 Int_t jrow=irow;
928 if (isec>=36) jrow+=63;
af6a50bb 929 for (UInt_t ipad=0;ipad<proc->GetNPads(isec,irow);ipad+=1){
930 proc->GetPositionGlobal(isec,irow,ipad, globalPos);
931 Double_t phi=TMath::ATan2(globalPos[1],globalPos[0]);
932 //if (phi<0) phi+=TMath::Pi()*2;
933 Int_t bin=his->FindBin(phi,jrow);
934 Float_t value= his->GetBinContent(bin);
935 calRoc->SetValue(irow,ipad,value);
936 }
937 }
938 }
939 return calPad;
940}
0e00d630 941
e52a66af 942AliTPCCalPad *AliTPCCalPad::MakePadFromTree(TTree * treePad, const char *query, const char* name, Bool_t doFast){
0e00d630 943 //
944 // make cal pad from the tree
945 //
99d61d46 946 if (!treePad){
947 ::Error("AliTPCCalPad::MakePadFromTree(TTree * treePad, const char *query, const char* name)","Input tree is missing");
948 return 0;
949 }
0e00d630 950 if (treePad->GetEntries()!=kNsec) return 0;
951 AliTPCCalPad * calPad= new AliTPCCalPad(name,name);
952 if (name) calPad->SetName(name);
e52a66af 953 if (!doFast){
954 for (Int_t iSec=0; iSec<72; iSec++){
955 AliTPCCalROC* calROC = calPad->GetCalROC(iSec);
956 UInt_t nchannels = (UInt_t)treePad->Draw(query,"1","goff",1,iSec);
957 if (nchannels!=calROC->GetNchannels()) {
958 ::Error("AliTPCCalPad::MakePad",TString::Format("%s\t:Wrong query sector\t%d\t%d",treePad->GetName(),iSec,nchannels).Data());
959 break;
960 }
961 for (UInt_t index=0; index<nchannels; index++) calROC->SetValue(index,treePad->GetV1()[index]);
962 }
963 }else{
964 UInt_t nchannelsTree = (UInt_t)treePad->Draw(query,"1","goff");
965 UInt_t nchannelsAll=0;
966 for (Int_t iSec=0; iSec<72; iSec++){
967 AliTPCCalROC* calROC = calPad->GetCalROC(iSec);
968 UInt_t nchannels=calROC->GetNchannels();
969 for (UInt_t index=0; index<nchannels; index++) {
970 if (nchannelsAll<=nchannelsTree)calROC->SetValue(index,treePad->GetV1()[nchannelsAll]);
971 nchannelsAll++;
972 }
973 }
974 if (nchannelsAll>nchannelsTree){
975 ::Error("AliTPCCalPad::MakePad",TString::Format("%s\t:Wrong query: cout mismatch\t%d\t%d",query, nchannelsAll,nchannelsTree).Data());
0e00d630 976 }
0e00d630 977 }
978 return calPad;
979}
980
981void AliTPCCalPad::AddFriend(TTree * treePad, const char *friendName, const char *fname){
982 //
983 //
984 //
985 TObjArray *fArray = new TObjArray(1);
986 fArray->AddLast(this);
987 this->SetName(friendName);
988 AliTPCCalibViewer::MakeTree(fname, fArray,0);
989 TFile * f = TFile::Open(fname);
990 TTree * tree = (TTree*)f->Get("calPads");
991 treePad->AddFriend(tree,friendName);
992 // tree->AddFriend(TString::Format("%s = calPads",friendName).Data(),fname);
993}