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