1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // TPC calibration class for parameters which are saved per pad //
21 // Each AliTPCCalPad consists of 72 AliTPCCalROC-objects //
23 ///////////////////////////////////////////////////////////////////////////////
25 #include "AliTPCCalPad.h"
26 #include "AliTPCCalROC.h"
27 #include <TObjArray.h>
32 #include "TTreeStream.h"
37 ClassImp(AliTPCCalPad)
39 //_____________________________________________________________________________
40 AliTPCCalPad::AliTPCCalPad():TNamed()
43 // AliTPCCalPad default constructor
46 for (Int_t isec = 0; isec < kNsec; isec++) {
52 //_____________________________________________________________________________
53 AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title)
57 // AliTPCCalPad constructor
59 for (Int_t isec = 0; isec < kNsec; isec++) {
60 fROC[isec] = new AliTPCCalROC(isec);
65 //_____________________________________________________________________________
66 AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c)
69 // AliTPCCalPad copy constructor
72 for (Int_t isec = 0; isec < kNsec; isec++) {
75 fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
79 //_____________________________________________________________________________
80 AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed()
83 // AliTPCCalPad default constructor
86 for (Int_t isec = 0; isec < kNsec; isec++) {
87 fROC[isec] = (AliTPCCalROC *)array->At(isec);
93 ///_____________________________________________________________________________
94 AliTPCCalPad::~AliTPCCalPad()
97 // AliTPCCalPad destructor
100 for (Int_t isec = 0; isec < kNsec; isec++) {
109 //_____________________________________________________________________________
110 AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c)
113 // Assignment operator
116 if (this != &c) ((AliTPCCalPad &) c).Copy(*this);
121 //_____________________________________________________________________________
122 void AliTPCCalPad::Copy(TObject &c) const
128 for (Int_t isec = 0; isec < kNsec; isec++) {
130 fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]);
137 void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){
139 // Set AliTPCCalROC copies values from 'roc'
140 // if sector == -1 the sector specified in 'roc' is used
141 // else sector specified in 'roc' is ignored and specified sector is filled
143 if (sector == -1) sector = roc->GetSector();
144 if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector);
145 for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++)
146 fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
151 //_____________________________________________________________________________
152 void AliTPCCalPad::Add(Float_t c1)
155 // add constant c1 to all channels of all ROCs
158 for (Int_t isec = 0; isec < kNsec; isec++) {
165 //_____________________________________________________________________________
166 void AliTPCCalPad::Multiply(Float_t c1)
169 // multiply each channel of all ROCs with c1
171 for (Int_t isec = 0; isec < kNsec; isec++) {
173 fROC[isec]->Multiply(c1);
178 //_____________________________________________________________________________
179 void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
182 // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
185 for (Int_t isec = 0; isec < kNsec; isec++) {
187 fROC[isec]->Add(pad->GetCalROC(isec),c1);
192 //_____________________________________________________________________________
193 void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
196 // multiply each channel of all ROCs with the coresponding channel of 'pad'
199 for (Int_t isec = 0; isec < kNsec; isec++) {
201 fROC[isec]->Multiply(pad->GetCalROC(isec));
206 //_____________________________________________________________________________
207 void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
210 // divide each channel of all ROCs by the coresponding channel of 'pad'
213 for (Int_t isec = 0; isec < kNsec; isec++) {
215 fROC[isec]->Divide(pad->GetCalROC(isec));
220 //_____________________________________________________________________________
221 TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
228 for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
229 TGraph * graph = new TGraph(npoints);
231 for (Int_t isec=0;isec<72;isec++){
232 if (!fROC[isec]) continue;
233 if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
234 if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
235 if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
239 graph->GetXaxis()->SetTitle("Sector");
241 graph->GetYaxis()->SetTitle("Mean");
242 graph->SetMarkerStyle(22);
245 graph->GetYaxis()->SetTitle("Median");
246 graph->SetMarkerStyle(22);
249 graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
250 graph->SetMarkerStyle(24);
256 //_____________________________________________________________________________
257 Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
260 // Calculates mean and RMS of all ROCs
262 Double_t sum = 0, sum2 = 0, n=0, val=0;
263 for (Int_t isec = 0; isec < kNsec; isec++) {
264 AliTPCCalROC *calRoc = fROC[isec];
266 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
267 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
268 val = calRoc->GetValue(irow,ipad);
278 Double_t mean = sum*n1;
279 rms = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
284 //_____________________________________________________________________________
285 Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
288 // return mean of the mean of all ROCs
292 for (Int_t isec = 0; isec < kNsec; isec++) {
293 AliTPCCalROC *calRoc = fROC[isec];
295 AliTPCCalROC* outlierROC = 0;
296 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
297 arr[n] = calRoc->GetMean(outlierROC);
301 return TMath::Mean(n,arr);
304 //_____________________________________________________________________________
305 Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
308 // return mean of the RMS of all ROCs
312 for (Int_t isec = 0; isec < kNsec; isec++) {
313 AliTPCCalROC *calRoc = fROC[isec];
315 AliTPCCalROC* outlierROC = 0;
316 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
317 arr[n] = calRoc->GetRMS(outlierROC);
321 return TMath::Mean(n,arr);
324 //_____________________________________________________________________________
325 Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
328 // return mean of the median of all ROCs
332 for (Int_t isec = 0; isec < kNsec; isec++) {
333 AliTPCCalROC *calRoc = fROC[isec];
335 AliTPCCalROC* outlierROC = 0;
336 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
337 arr[n] = calRoc->GetMedian(outlierROC);
341 return TMath::Mean(n,arr);
344 //_____________________________________________________________________________
345 Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
348 // return mean of the LTM and sigma of all ROCs
350 Double_t arrm[kNsec];
351 Double_t arrs[kNsec];
355 for (Int_t isec = 0; isec < kNsec; isec++) {
356 AliTPCCalROC *calRoc = fROC[isec];
358 if ( sigma ) sTemp=arrs+n;
359 AliTPCCalROC* outlierROC = 0;
360 if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
361 arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
365 if ( sigma ) *sigma = TMath::Mean(n,arrs);
366 return TMath::Mean(n,arrm);
369 //_____________________________________________________________________________
370 TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){
373 // type -1 = user defined range
374 // 0 = nsigma cut nsigma=min
379 Float_t mean = GetMean();
380 Float_t sigma = GetRMS();
381 Float_t nsigma = TMath::Abs(min);
382 min = mean-nsigma*sigma;
383 max = mean+nsigma*sigma;
387 Float_t mean = GetMedian();
394 // LTM mean +- nsigma
397 Float_t mean = GetLTM(&sigma,max);
404 sprintf(name,"%s Pad 1D",GetTitle());
405 TH1F * his = new TH1F(name,name,100, min,max);
406 for (Int_t isec = 0; isec < kNsec; isec++) {
408 for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
409 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
410 for (UInt_t ipad=0; ipad<npads; ipad++){
411 his->Fill(fROC[isec]->GetValue(irow,ipad));
419 //_____________________________________________________________________________
420 TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
423 // side - specify the side A = 0 C = 1
424 // type - used types of determination of boundaries in z
426 Float_t kEpsilon = 0.000000000001;
427 TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
428 AliTPCROC * roc = AliTPCROC::Instance();
429 for (Int_t isec=0; isec<72; isec++){
430 if (side==0 && isec%36>=18) continue;
431 if (side>0 && isec%36<18) continue;
433 AliTPCCalROC * calRoc = fROC[isec];
434 for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
435 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
436 if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
438 roc->GetPositionGlobal(isec,irow,ipad,xyz);
439 Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
440 Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
441 Float_t value = calRoc->GetValue(irow,ipad);
442 his->SetBinContent(binx,biny,value);
446 his->SetXTitle("x (cm)");
447 his->SetYTitle("y (cm)");
452 AliTPCCalPad* 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 {
454 // Loops over all AliTPCCalROCs and performs a localFit in each ROC
455 // AliTPCCalPad with fit-data is returned
456 // rowRadius and padRadius specifies a window around a given pad in one sector.
457 // The data of this window are fitted with a parabolic function.
458 // This function is evaluated at the pad's position.
459 // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
460 // rowRadius - radius - rows to be used for smoothing
461 // padradius - radius - pads to be used for smoothing
462 // ROCoutlier - map of outliers - pads not to be used for local smoothing
463 // robust - robust method of fitting - (much slower)
464 // chi2Threshold: Threshold for chi2 when EvalRobust is called
465 // robustFraction: Fraction of data that will be used in EvalRobust
468 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
469 for (Int_t isec = 0; isec < 72; isec++){
470 if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
472 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
474 pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
480 AliTPCCalPad* 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){
482 // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
483 // AliTPCCalPad with fit-data is returned
484 // chi2Threshold: Threshold for chi2 when EvalRobust is called
485 // robustFraction: Fraction of data that will be used in EvalRobust
486 // chi2Threshold: Threshold for chi2 when EvalRobust is called
487 // robustFraction: Fraction of data that will be used in EvalRobust
488 // err: error of the data points
489 // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
491 AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
492 TVectorD fitParam(0);
493 TMatrixD covMatrix(0,0);
495 for (Int_t isec = 0; isec < 72; isec++){
497 GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
499 GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
501 AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec);
504 if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec);
505 if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec);
511 void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, Int_t fitType, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){
513 // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side
514 // fitType == 0: fit plane function
515 // fitType == 1: fit parabolic function
516 // PadOutliers - pads with value !=0 are not used in fitting procedure
517 // chi2Threshold: Threshold for chi2 when EvalRobust is called
518 // robustFraction: Fraction of data that will be used in EvalRobust
520 TLinearFitter* fitterGA = 0;
521 TLinearFitter* fitterGC = 0;
524 fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
525 fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
528 fitterGA = new TLinearFitter(3,"x0++x1++x2");
529 fitterGC = new TLinearFitter(3,"x0++x1++x2");
531 fitterGA->StoreData(kTRUE);
532 fitterGC->StoreData(kTRUE);
533 fitterGA->ClearPoints();
534 fitterGC->ClearPoints();
539 Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position
540 Float_t lx, ly; // pads position
542 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
544 // loop over all sectors and pads and read data into fitterGA and fitterGC
547 fitParamSideA.ResizeTo(6);
548 fitParamSideC.ResizeTo(6);
549 covMatrixSideA.ResizeTo(6,6);
550 covMatrixSideC.ResizeTo(6,6);
551 for (UInt_t isec = 0; isec<72; isec++){
552 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
553 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
555 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
564 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
565 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
566 // sector 0 - 17: IROC, A
567 // sector 18 - 35: IROC, C
568 // sector 36 - 53: OROC, A
569 // sector 54 - 71: CROC, C
570 if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A
572 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
576 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
585 fitParamSideA.ResizeTo(3);
586 fitParamSideC.ResizeTo(3);
587 covMatrixSideA.ResizeTo(3,3);
588 covMatrixSideC.ResizeTo(3,3);
590 for (UInt_t isec = 0; isec<72; isec++){
591 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
592 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
594 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
600 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
601 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
602 // sector 0 - 17: IROC, A
603 // sector 18 - 35: IROC, C
604 // sector 36 - 53: OROC, A
605 // sector 54 - 71: CROC, C
606 if (isec <= 17 || (isec >= 36 && isec <= 53)) {
609 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
614 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
624 fitterGA->GetParameters(fitParamSideA);
625 fitterGC->GetParameters(fitParamSideC);
626 fitterGA->GetCovarianceMatrix(covMatrixSideA);
627 fitterGC->GetCovarianceMatrix(covMatrixSideC);
629 chi2SideA = fitterGA->GetChisquare()/(npointsA-6.);
630 chi2SideC = fitterGC->GetChisquare()/(npointsC-6.);
633 chi2SideA = fitterGA->GetChisquare()/(npointsA-3.);
634 chi2SideC = fitterGC->GetChisquare()/(npointsC-3.);
636 if (robust && chi2SideA > chi2Threshold) {
637 // std::cout << "robust fitter called... " << std::endl;
638 fitterGA->EvalRobust(robustFraction);
639 fitterGA->GetParameters(fitParamSideA);
641 if (robust && chi2SideC > chi2Threshold) {
642 // std::cout << "robust fitter called... " << std::endl;
643 fitterGC->EvalRobust(robustFraction);
644 fitterGC->GetParameters(fitParamSideC);