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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Class providing the calculation of derived quantities (mean,rms,fits,...) //
20 // of calibration entries //
25 ////////////////////////////////////////////////////////////////////////////////
29 #include <TObjArray.h>
32 #include <TDirectory.h>
34 #include <TGraphErrors.h>
35 #include <AliCDBStorage.h>
36 #include <AliDCSSensorArray.h>
37 #include <AliTPCSensorTempArray.h>
38 #include <AliDCSSensor.h>
40 #include <AliCDBEntry.h>
41 #include <AliCDBManager.h>
43 #include <AliSplineFit.h>
44 #include "AliTPCcalibDB.h"
45 #include "AliTPCCalPad.h"
46 #include "AliTPCCalROC.h"
47 #include "AliTPCROC.h"
48 #include "AliTPCmapper.h"
49 #include "AliTPCParam.h"
50 #include "AliTPCCalibRaw.h"
51 #include "AliTPCPreprocessorOnline.h"
52 #include "AliTPCdataQA.h"
54 #include "AliTPCcalibDButil.h"
55 #include "AliTPCCalibVdrift.h"
56 #include "AliMathBase.h"
57 #include "AliRelAlignerKalman.h"
62 const Float_t kAlmost0=1.e-30;
64 ClassImp(AliTPCcalibDButil)
65 AliTPCcalibDButil::AliTPCcalibDButil() :
73 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
85 fRefPedestalMasked(0x0),
89 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
90 fRefPulserMasked(0x0),
97 fRefALTROAcqStart(0x0),
98 fRefALTROAcqStop(0x0),
103 fMapper(new AliTPCmapper(0x0)),
104 fNpulserOutliers(-1),
106 fCETmaxLimitAbs(1.5),
107 fPulTmaxLimitAbs(1.5),
110 fRuns(0), // run list with OCDB info
111 fRunsStart(0), // start time for given run
112 fRunsStop(0) // stop time for given run
118 //_____________________________________________________________________________________
119 AliTPCcalibDButil::~AliTPCcalibDButil()
124 delete fPulserOutlier;
125 delete fRefPulserOutlier;
127 if (fRefPadNoise) delete fRefPadNoise;
128 if (fRefPedestals) delete fRefPedestals;
129 if (fRefPedestalMasked) delete fRefPedestalMasked;
130 if (fRefPulserTmean) delete fRefPulserTmean;
131 if (fRefPulserTrms) delete fRefPulserTrms;
132 if (fRefPulserQmean) delete fRefPulserQmean;
133 if (fRefPulserMasked) delete fRefPulserMasked;
134 if (fRefCETmean) delete fRefCETmean;
135 if (fRefCETrms) delete fRefCETrms;
136 if (fRefCEQmean) delete fRefCEQmean;
137 if (fRefCEMasked) delete fRefCEMasked;
138 if (fRefALTROFPED) delete fRefALTROFPED;
139 if (fRefALTROZsThr) delete fRefALTROZsThr;
140 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
141 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
142 if (fRefALTROMasked) delete fRefALTROMasked;
143 if (fRefCalibRaw) delete fRefCalibRaw;
144 if (fCurrentRefMap) delete fCurrentRefMap;
146 //_____________________________________________________________________________________
147 void AliTPCcalibDButil::UpdateFromCalibDB()
150 // Update pointers from calibDB
152 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
153 fCalibDB->UpdateNonRec(); // load all infromation now
154 fPadNoise=fCalibDB->GetPadNoise();
155 fPedestals=fCalibDB->GetPedestals();
156 fPulserTmean=fCalibDB->GetPulserTmean();
157 fPulserTrms=fCalibDB->GetPulserTrms();
158 fPulserQmean=fCalibDB->GetPulserQmean();
159 fCETmean=fCalibDB->GetCETmean();
160 fCETrms=fCalibDB->GetCETrms();
161 fCEQmean=fCalibDB->GetCEQmean();
162 fALTROMasked=fCalibDB->GetALTROMasked();
163 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
164 fCalibRaw=fCalibDB->GetCalibRaw();
165 fDataQA=fCalibDB->GetDataQA();
166 UpdatePulserOutlierMap();
167 // SetReferenceRun();
168 // UpdateRefDataFromOCDB();
170 //_____________________________________________________________________________________
171 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
172 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
175 // Process the CE data for this run
176 // the return TVectorD arrays contian the results of the fit
177 // noutliersCE contains the number of pads marked as outliers,
178 // not including masked and edge pads
181 //retrieve CE and ALTRO data
183 TString fitString(fitFormula);
184 fitString.ReplaceAll("++","#");
185 Int_t ndim=fitString.CountChar('#')+2;
186 fitResultsA.ResizeTo(ndim);
187 fitResultsC.ResizeTo(ndim);
196 if (outCE) out=outCE;
197 else out=new AliTPCCalPad("outCE","outCE");
198 AliTPCCalROC *rocMasked=0x0;
199 //loop over all channels
200 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
201 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
202 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
203 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
205 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
209 //add time offset to IROCs
210 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
211 rocData->Add(fIrocTimeOffset);
213 UInt_t nrows=rocData->GetNrows();
214 for (UInt_t irow=0;irow<nrows;++irow){
215 UInt_t npads=rocData->GetNPads(irow);
216 for (UInt_t ipad=0;ipad<npads;++ipad){
217 rocOut->SetValue(irow,ipad,0);
218 //exclude masked pads
219 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
220 rocOut->SetValue(irow,ipad,1);
223 //exclude first two rows in IROC and last two rows in OROC
225 if (irow<2) rocOut->SetValue(irow,ipad,1);
227 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
230 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
231 Float_t valTmean=rocData->GetValue(irow,ipad);
232 //exclude values that are exactly 0
233 if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
234 rocOut->SetValue(irow,ipad,1);
237 // exclude channels with too large variations
238 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
239 rocOut->SetValue(irow,ipad,1);
247 Float_t chi2Af,chi2Cf;
248 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
251 if (!outCE) delete out;
253 //_____________________________________________________________________________________
254 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
255 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
256 Float_t &driftTimeA, Float_t &driftTimeC )
259 // Calculate statistical information from the CE graphs for drift time and charge
263 vecTEntries.ResizeTo(72);
264 vecTMean.ResizeTo(72);
265 vecTRMS.ResizeTo(72);
266 vecTMedian.ResizeTo(72);
267 vecQEntries.ResizeTo(72);
268 vecQMean.ResizeTo(72);
269 vecQRMS.ResizeTo(72);
270 vecQMedian.ResizeTo(72);
281 TObjArray *arrT=fCalibDB->GetCErocTtime();
282 TObjArray *arrQ=fCalibDB->GetCErocQtime();
284 for (Int_t isec=0;isec<74;++isec){
285 TGraph *gr=(TGraph*)arrT->At(isec);
288 Int_t npoints = gr->GetN();
289 values.ResizeTo(npoints);
291 //skip first points, theres always some problems with finding the CE position
292 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
293 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
294 values[nused]=gr->GetY()[ipoint];
299 if (isec<72) vecTEntries[isec]= nused;
302 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
303 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
304 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
305 } else if (isec==72){
306 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
307 } else if (isec==73){
308 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
314 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
315 TGraph *gr=(TGraph*)arrQ->At(isec);
318 Int_t npoints = gr->GetN();
319 values.ResizeTo(npoints);
321 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
322 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
323 values[nused]=gr->GetY()[ipoint];
328 vecQEntries[isec]= nused;
330 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
331 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
332 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
338 //_____________________________________________________________________________________
339 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
340 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
341 Int_t &nonMaskedZero, Int_t &nNaN)
344 // process noise data
345 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
346 // OROCs small pads [2] and OROCs large pads [3]
347 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
348 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
351 //set proper size and reset
352 const UInt_t infoSize=4;
353 vNoiseMean.ResizeTo(infoSize);
354 vNoiseMeanSenRegions.ResizeTo(infoSize);
355 vNoiseRMS.ResizeTo(infoSize);
356 vNoiseRMSSenRegions.ResizeTo(infoSize);
358 vNoiseMeanSenRegions.Zero();
360 vNoiseRMSSenRegions.Zero();
364 TVectorD c(infoSize);
365 TVectorD cs(infoSize);
369 //retrieve noise and ALTRO data
370 if (!fPadNoise) return;
371 AliTPCCalROC *rocMasked=0x0;
372 //create IROC, OROC1, OROC2 and sensitive region masks
373 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
374 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
375 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
376 UInt_t nrows=noiseROC->GetNrows();
377 for (UInt_t irow=0;irow<nrows;++irow){
378 UInt_t npads=noiseROC->GetNPads(irow);
379 for (UInt_t ipad=0;ipad<npads;++ipad){
380 //don't use masked channels;
381 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
382 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
384 if (noiseVal<kAlmost0) {
389 if ( !(noiseVal<10000000) ){
390 AliInfo(Form("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal));
394 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
395 Int_t masksen=1; // sensitive pards are not masked (0)
396 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
397 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
399 if (irow>19&&irow<46){
400 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
403 vNoiseMean[type]+=noiseVal;
404 vNoiseRMS[type]+=noiseVal*noiseVal;
407 vNoiseMeanSenRegions[type]+=noiseVal;
408 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
413 //define sensive regions
414 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
416 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
417 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
419 if ((Int_t)irow<par.GetNRowUp1()){
422 vNoiseMean[type]+=noiseVal;
423 vNoiseRMS[type]+=noiseVal*noiseVal;
426 vNoiseMeanSenRegions[type]+=noiseVal;
427 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
433 vNoiseMean[type]+=noiseVal;
434 vNoiseRMS[type]+=noiseVal*noiseVal;
437 vNoiseMeanSenRegions[type]+=noiseVal;
438 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
445 vNoiseMean[type]+=noiseVal;
446 vNoiseRMS[type]+=noiseVal*noiseVal;
449 vNoiseMeanSenRegions[type]+=noiseVal;
450 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
455 }//end loop sectors (rocs)
457 //calculate mean and RMS
458 const Double_t verySmall=0.0000000001;
459 for (UInt_t i=0;i<infoSize;++i){
466 AliInfo(Form("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]));
467 mean=vNoiseMean[i]/c[i];
469 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
474 if (cs[i]>verySmall){
475 meanSen=vNoiseMeanSenRegions[i]/cs[i];
476 rmsSen=vNoiseRMSSenRegions[i];
477 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
479 vNoiseMeanSenRegions[i]=meanSen;
480 vNoiseRMSSenRegions[i]=rmsSen;
484 //_____________________________________________________________________________________
485 void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
491 // vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
495 const UInt_t infoSize = 72;
496 //reset counters to error number
497 vQaOcc.ResizeTo(infoSize);
499 vQaQtot.ResizeTo(infoSize);
501 vQaQmax.ResizeTo(infoSize);
504 //retrieve pulser and ALTRO data
508 AliInfo("No QA data");
511 if (fDataQA->GetEventCounter()<=0) {
513 AliInfo("No QA data");
514 return; // no data processed
519 TVectorD normOcc(infoSize);
520 TVectorD normQ(infoSize);
522 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
524 AliInfo(Form("Sector %d\n", isec));
525 AliTPCCalROC* occupancyROC = fDataQA->GetNoThreshold()->GetCalROC(isec);
526 AliTPCCalROC* nclusterROC = fDataQA->GetNLocalMaxima()->GetCalROC(isec);
527 AliTPCCalROC* qROC = fDataQA->GetMeanCharge()->GetCalROC(isec);
528 AliTPCCalROC* qmaxROC = fDataQA->GetMaxCharge()->GetCalROC(isec);
529 if (!occupancyROC) continue;
530 if (!nclusterROC) continue;
532 if (!qmaxROC) continue;
534 const UInt_t nchannels=occupancyROC->GetNchannels();
536 AliInfo(Form("Nchannels %d\n", nchannels));
538 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
540 vQaOcc[isec] += occupancyROC->GetValue(ichannel);
543 Float_t nClusters = nclusterROC->GetValue(ichannel);
544 normQ[isec] += nClusters;
545 vQaQtot[isec]+=nClusters*qROC->GetValue(ichannel);
546 vQaQmax[isec]+=nClusters*qmaxROC->GetValue(ichannel);
550 //calculate mean values
551 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
553 if (normOcc[isec]>0) vQaOcc[isec] /= normOcc[isec];
554 else vQaOcc[isec] = 0;
557 vQaQtot[isec] /= normQ[isec];
558 vQaQmax[isec] /= normQ[isec];
567 //_____________________________________________________________________________________
568 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
571 // Process the Pulser information
572 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
575 const UInt_t infoSize=4;
576 //reset counters to error number
577 vMeanTime.ResizeTo(infoSize);
580 TVectorD c(infoSize);
581 //retrieve pulser and ALTRO data
582 if (!fPulserTmean) return;
585 AliTPCCalROC *rocOut=0x0;
586 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
587 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
588 if (!tmeanROC) continue;
589 rocOut=fPulserOutlier->GetCalROC(isec);
590 UInt_t nchannels=tmeanROC->GetNchannels();
591 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
592 if (rocOut && rocOut->GetValue(ichannel)) continue;
593 Float_t val=tmeanROC->GetValue(ichannel);
595 vMeanTime[type]+=val;
600 for (UInt_t itype=0; itype<infoSize; ++itype){
601 if (c[itype]>0) vMeanTime[itype]/=c[itype];
602 else vMeanTime[itype]=0;
605 //_____________________________________________________________________________________
606 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
609 // Get Values from ALTRO configuration data
612 if (!fALTROMasked) return;
614 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
615 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
616 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
617 if (rocMasked->GetValue(ichannel)) ++nMasked;
621 //_____________________________________________________________________________________
622 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
625 // Proces Goofie values, return statistical information of the currently set goofieArray
626 // The meaning of the entries are given below
628 1 TPC_ANODE_I_A00_STAT
630 3 TPC_DVM_DriftVelocity
635 8 TPC_DVM_NumberOfSparks
636 9 TPC_DVM_PeakAreaFar
637 10 TPC_DVM_PeakAreaNear
638 11 TPC_DVM_PeakPosFar
639 12 TPC_DVM_PeakPosNear
645 18 TPC_DVM_TemperatureS1
649 vecEntries.ResizeTo(nsensors);
650 vecMedian.ResizeTo(nsensors);
651 vecMean.ResizeTo(nsensors);
652 vecRMS.ResizeTo(nsensors);
659 Double_t kEpsilon=0.0000000001;
660 Double_t kBig=100000000000.;
661 Int_t nsensors = fGoofieArray->NumSensors();
662 vecEntries.ResizeTo(nsensors);
663 vecMedian.ResizeTo(nsensors);
664 vecMean.ResizeTo(nsensors);
665 vecRMS.ResizeTo(nsensors);
667 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
668 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
669 if (gsensor && gsensor->GetGraph()){
670 Int_t npoints = gsensor->GetGraph()->GetN();
672 values.ResizeTo(npoints);
674 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
675 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
676 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
677 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
682 vecEntries[isensor]= nused;
684 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
685 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
686 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
691 //_____________________________________________________________________________________
692 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
695 // check the variations of the pedestal data to the reference pedestal data
696 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
699 TVectorF vThres(npar); //thresholds
700 Int_t nActive=0; //number of active channels
702 //reset and set thresholds
703 pedestalDeviations.ResizeTo(npar);
704 for (Int_t i=0;i<npar;++i){
705 pedestalDeviations.GetMatrixArray()[i]=0;
706 vThres.GetMatrixArray()[i]=(i+1)*.5;
708 //check all needed data is available
709 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
710 //loop over all channels
711 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
712 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
713 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
714 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
715 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
716 UInt_t nrows=mROC->GetNrows();
717 for (UInt_t irow=0;irow<nrows;++irow){
718 UInt_t npads=mROC->GetNPads(irow);
719 for (UInt_t ipad=0;ipad<npads;++ipad){
720 //don't use masked channels;
721 if (mROC ->GetValue(irow,ipad)) continue;
722 if (mRefROC->GetValue(irow,ipad)) continue;
723 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
724 for (Int_t i=0;i<npar;++i){
725 if (deviation>vThres[i])
726 ++pedestalDeviations.GetMatrixArray()[i];
733 for (Int_t i=0;i<npar;++i){
734 pedestalDeviations.GetMatrixArray()[i]/=nActive;
738 //_____________________________________________________________________________________
739 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
742 // check the variations of the noise data to the reference noise data
743 // thresholds are 5, 10, 15 and 20 percent respectively.
746 TVectorF vThres(npar); //thresholds
747 Int_t nActive=0; //number of active channels
749 //reset and set thresholds
750 noiseDeviations.ResizeTo(npar);
751 for (Int_t i=0;i<npar;++i){
752 noiseDeviations.GetMatrixArray()[i]=0;
753 vThres.GetMatrixArray()[i]=(i+1)*.05;
755 //check all needed data is available
756 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
757 //loop over all channels
758 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
759 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
760 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
761 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
762 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
763 UInt_t nrows=mROC->GetNrows();
764 for (UInt_t irow=0;irow<nrows;++irow){
765 UInt_t npads=mROC->GetNPads(irow);
766 for (UInt_t ipad=0;ipad<npads;++ipad){
767 //don't use masked channels;
768 if (mROC ->GetValue(irow,ipad)) continue;
769 if (mRefROC->GetValue(irow,ipad)) continue;
770 if (nRefROC->GetValue(irow,ipad)==0) continue;
771 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
772 for (Int_t i=0;i<npar;++i){
773 if (deviation>vThres[i])
774 ++noiseDeviations.GetMatrixArray()[i];
781 for (Int_t i=0;i<npar;++i){
782 noiseDeviations.GetMatrixArray()[i]/=nActive;
786 //_____________________________________________________________________________________
787 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
788 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
791 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
792 // thresholds are .5, 1, 5 and 10 percent respectively.
796 TVectorF vThres(npar); //thresholds
797 Int_t nActive=0; //number of active channels
799 //reset and set thresholds
800 pulserQdeviations.ResizeTo(npar);
801 for (Int_t i=0;i<npar;++i){
802 pulserQdeviations.GetMatrixArray()[i]=0;
807 vThres.GetMatrixArray()[0]=.005;
808 vThres.GetMatrixArray()[1]=.01;
809 vThres.GetMatrixArray()[2]=.05;
810 vThres.GetMatrixArray()[3]=.1;
811 //check all needed data is available
812 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
814 UpdateRefPulserOutlierMap();
815 //loop over all channels
816 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
817 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
818 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
819 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
820 // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
821 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
822 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
823 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
824 Float_t ptmean=ptROC->GetMean(oROC);
825 UInt_t nrows=mROC->GetNrows();
826 for (UInt_t irow=0;irow<nrows;++irow){
827 UInt_t npads=mROC->GetNPads(irow);
828 for (UInt_t ipad=0;ipad<npads;++ipad){
829 //don't use masked channels;
830 if (mROC ->GetValue(irow,ipad)) continue;
831 if (mRefROC->GetValue(irow,ipad)) continue;
832 //don't user edge pads
833 if (ipad==0||ipad==npads-1) continue;
835 Float_t pq=pqROC->GetValue(irow,ipad);
836 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
837 Float_t pt=ptROC->GetValue(irow,ipad);
838 // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
840 Float_t deviation=TMath::Abs(pqRef)>1e-20?TMath::Abs(pq/pqRef-1):-999;
841 for (Int_t i=0;i<npar;++i){
842 if (deviation>vThres[i])
843 ++pulserQdeviations.GetMatrixArray()[i];
845 if (pqRef>11&&pq<11) ++npadsOffAdd;
848 if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
854 for (Int_t i=0;i<npar;++i){
855 pulserQdeviations.GetMatrixArray()[i]/=nActive;
860 //_____________________________________________________________________________________
861 void AliTPCcalibDButil::UpdatePulserOutlierMap()
864 // Update the outlier map of the pulser data
866 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
868 //_____________________________________________________________________________________
869 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
872 // Update the outlier map of the pulser reference data
874 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
876 //_____________________________________________________________________________________
877 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
880 // Create a map that contains outliers from the Pulser calibration data.
881 // The outliers include masked channels, edge pads and pads with
882 // too large timing and charge variations.
883 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
884 // those do not contain masked and edge pads
888 pulOut->Multiply(0.);
892 AliTPCCalROC *rocMasked=0x0;
896 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
897 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
898 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
899 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
900 if (!tmeanROC||!qmeanROC) {
901 //reset outliers in this ROC
902 outROC->Multiply(0.);
905 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
907 // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
908 // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
909 UInt_t nrows=tmeanROC->GetNrows();
910 for (UInt_t irow=0;irow<nrows;++irow){
911 UInt_t npads=tmeanROC->GetNPads(irow);
912 for (UInt_t ipad=0;ipad<npads;++ipad){
913 Int_t outlier=0,masked=0;
914 Float_t q=qmeanROC->GetValue(irow,ipad);
915 Float_t t=tmeanROC->GetValue(irow,ipad);
916 //masked channels are outliers
917 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
918 //edge pads are outliers
919 if (ipad==0||ipad==npads-1) masked=1;
920 //channels with too large charge or timing deviation from the meadian are outliers
921 // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
922 if (q<fPulQminLimit && !masked) outlier=1;
924 if ( !(q<10000000) || !(t<10000000)) outlier=1;
925 outROC->SetValue(irow,ipad,outlier+masked);
926 fNpulserOutliers+=outlier;
931 //_____________________________________________________________________________________
932 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
935 // Create pad time0 object from pulser and/or CE data, depending on the selected model
936 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
937 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
938 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
940 // In case model 2 is invoked - gy arival time gradient is also returned
944 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
945 // decide between different models
946 if (model==0||model==1){
948 if (model==1) ProcessPulser(vMean);
949 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
950 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
951 if (!rocPulTmean) continue;
952 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
953 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
954 Float_t mean=rocPulTmean->GetMean(rocOut);
955 //treat case where a whole partition is masked
956 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
961 UInt_t nrows=rocTime0->GetNrows();
962 for (UInt_t irow=0;irow<nrows;++irow){
963 UInt_t npads=rocTime0->GetNPads(irow);
964 for (UInt_t ipad=0;ipad<npads;++ipad){
965 Float_t time=rocPulTmean->GetValue(irow,ipad);
966 //in case of an outlier pad use the mean of the altro values.
967 //This should be the most precise guess in that case.
968 if (rocOut->GetValue(irow,ipad)) {
969 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
970 if ( TMath::Abs(time)<kAlmost0 ) time=mean;
972 Float_t val=time-mean;
973 rocTime0->SetValue(irow,ipad,val);
977 } else if (model==2){
978 Double_t pgya,pgyc,pchi2a,pchi2c;
979 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
980 fCETmean->Add(padPulser,-1.);
982 AliTPCCalPad outCE("outCE","outCE");
984 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
985 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
986 // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
987 if (!padFit) { delete padPulser; return 0;}
990 fCETmean->Add(padPulser,1.);
991 padTime0->Add(fCETmean);
992 padTime0->Add(padFit,-1);
997 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
998 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
999 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
1000 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
1001 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
1002 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
1003 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
1004 Float_t mean=rocPulTmean->GetMean(rocOutPul);
1005 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
1006 UInt_t nrows=rocTime0->GetNrows();
1007 for (UInt_t irow=0;irow<nrows;++irow){
1008 UInt_t npads=rocTime0->GetNPads(irow);
1009 for (UInt_t ipad=0;ipad<npads;++ipad){
1010 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
1011 if (rocOutCE->GetValue(irow,ipad)){
1012 Float_t valOut=rocCEfit->GetValue(irow,ipad);
1013 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
1014 rocTime0->SetValue(irow,ipad,valOut);
1022 Double_t median = padTime0->GetMedian();
1023 padTime0->Add(-median); // normalize to median
1026 //_____________________________________________________________________________________
1027 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
1030 // GetMeanAlto information
1032 if (roc==0) return 0.;
1033 const Int_t sector=roc->GetSector();
1034 AliTPCROC *tpcRoc=AliTPCROC::Instance();
1035 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
1039 //loop over a small range around the requested pad (+-10 rows/pads)
1040 for (Int_t irow=row-10;irow<row+10;++irow){
1041 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
1042 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
1043 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
1044 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
1045 if (altroRoc!=altroCurr) continue;
1046 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
1047 Float_t val=roc->GetValue(irow,ipad);
1055 //_____________________________________________________________________________________
1056 void AliTPCcalibDButil::SetRefFile(const char* filename)
1059 // load cal pad objects form the reference file
1061 TDirectory *currDir=gDirectory;
1063 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
1064 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
1066 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
1067 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
1068 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
1070 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
1071 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
1072 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
1074 // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
1075 // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
1076 // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
1077 // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
1078 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
1082 //_____________________________________________________________________________________
1083 void AliTPCcalibDButil::UpdateRefDataFromOCDB()
1086 // set reference data from OCDB Reference map
1089 AliWarning("Referenc map not set!");
1094 AliCDBEntry* entry = 0x0;
1095 Bool_t hasAnyChanged=kFALSE;
1098 cdbPath="TPC/Calib/Pedestals";
1099 if (HasRefChanged(cdbPath.Data())){
1100 hasAnyChanged=kTRUE;
1101 //delete old entries
1102 if (fRefPedestals) delete fRefPedestals;
1103 if (fRefPedestalMasked) delete fRefPedestalMasked;
1104 fRefPedestals=fRefPedestalMasked=0x0;
1106 entry=GetRefEntry(cdbPath.Data());
1108 entry->SetOwner(kTRUE);
1109 fRefPedestals=GetRefCalPad(entry);
1111 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1116 cdbPath="TPC/Calib/PadNoise";
1117 if (HasRefChanged(cdbPath.Data())){
1118 hasAnyChanged=kTRUE;
1120 if (fRefPadNoise) delete fRefPadNoise;
1123 entry=GetRefEntry(cdbPath.Data());
1125 entry->SetOwner(kTRUE);
1126 fRefPadNoise=GetRefCalPad(entry);
1132 cdbPath="TPC/Calib/Pulser";
1133 if (HasRefChanged(cdbPath.Data())){
1134 hasAnyChanged=kTRUE;
1135 //delete old entries
1136 if (fRefPulserTmean) delete fRefPulserTmean;
1137 if (fRefPulserTrms) delete fRefPulserTrms;
1138 if (fRefPulserQmean) delete fRefPulserQmean;
1139 if (fRefPulserMasked) delete fRefPulserMasked;
1140 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1142 entry=GetRefEntry(cdbPath.Data());
1144 entry->SetOwner(kTRUE);
1145 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1146 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1147 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1149 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1154 cdbPath="TPC/Calib/CE";
1155 if (HasRefChanged(cdbPath.Data())){
1156 hasAnyChanged=kTRUE;
1157 //delete old entries
1158 if (fRefCETmean) delete fRefCETmean;
1159 if (fRefCETrms) delete fRefCETrms;
1160 if (fRefCEQmean) delete fRefCEQmean;
1161 if (fRefCEMasked) delete fRefCEMasked;
1162 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1164 entry=GetRefEntry(cdbPath.Data());
1166 entry->SetOwner(kTRUE);
1167 fRefCETmean=GetRefCalPad(entry,"CETmean");
1168 fRefCETrms=GetRefCalPad(entry,"CETrms");
1169 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1171 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1176 cdbPath="TPC/Calib/AltroConfig";
1177 if (HasRefChanged(cdbPath.Data())){
1178 hasAnyChanged=kTRUE;
1179 //delete old entries
1180 if (fRefALTROFPED) delete fRefALTROFPED;
1181 if (fRefALTROZsThr) delete fRefALTROZsThr;
1182 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1183 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1184 if (fRefALTROMasked) delete fRefALTROMasked;
1185 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1187 entry=GetRefEntry(cdbPath.Data());
1189 entry->SetOwner(kTRUE);
1190 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1191 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1192 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1193 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1194 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1201 cdbPath="TPC/Calib/Raw";
1202 if (HasRefChanged(cdbPath.Data())){
1203 hasAnyChanged=kTRUE;
1205 if (fRefCalibRaw) delete fRefCalibRaw;
1207 entry=GetRefEntry(cdbPath.Data());
1209 entry->SetOwner(kTRUE);
1210 TObjArray *arr=(TObjArray*)entry->GetObject();
1212 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1214 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1221 cdbPath="TPC/Calib/QA";
1222 if (HasRefChanged(cdbPath.Data())){
1223 hasAnyChanged=kTRUE;
1225 if (fRefDataQA) delete fRefDataQA;
1227 entry=GetRefEntry(cdbPath.Data());
1229 entry->SetOwner(kTRUE);
1230 fRefDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1232 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1234 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
1241 //update current reference maps
1243 if (fCurrentRefMap) delete fCurrentRefMap;
1244 fCurrentRefMap=(TMap*)fRefMap->Clone();
1247 //_____________________________________________________________________________________
1248 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1251 // TObjArray object type case
1252 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1254 AliTPCCalPad *pad=0x0;
1255 TObjArray *arr=(TObjArray*)entry->GetObject();
1257 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1260 pad=(AliTPCCalPad*)arr->FindObject(objName);
1262 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
1265 return (AliTPCCalPad*)pad->Clone();
1267 //_____________________________________________________________________________________
1268 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1271 // AliTPCCalPad object type case
1272 // cast object to a calPad and store it in 'pad'
1274 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1276 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1279 pad=(AliTPCCalPad*)pad->Clone();
1282 //_____________________________________________________________________________________
1283 AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1286 // set altro masked channel map for 'cdbPath'
1288 AliTPCCalPad* pad=0x0;
1289 const Int_t run=GetReferenceRun(cdbPath);
1291 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1294 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1296 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1299 pad=GetRefCalPad(entry,"Masked");
1300 if (pad) pad->SetNameTitle(name,name);
1301 entry->SetOwner(kTRUE);
1305 //_____________________________________________________________________________________
1306 void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1308 // Get Reference map
1310 if (run<0) run=fCalibDB->GetRun();
1311 TString cdbPath="TPC/Calib/Ref";
1312 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1314 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
1318 entry->SetOwner(kTRUE);
1319 fRefMap=(TMap*)(entry->GetObject());
1320 AliCDBId &id=entry->GetId();
1321 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1323 //_____________________________________________________________________________________
1324 Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1327 // check whether a reference cdb entry has changed
1329 if (!fCurrentRefMap) return kTRUE;
1330 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1333 //_____________________________________________________________________________________
1334 AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1337 // get the reference AliCDBEntry for 'cdbPath'
1339 const Int_t run=GetReferenceRun(cdbPath);
1341 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1344 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1346 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1351 //_____________________________________________________________________________________
1352 Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
1354 // Get reference run number for the specified OCDB path
1356 if (!fCurrentRefMap) return -2;
1357 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1358 if (!str) return -2;
1359 return (Int_t)str->GetString().Atoi();
1361 //_____________________________________________________________________________________
1362 Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1364 // Get reference run number for the specified OCDB path
1366 if (!fRefMap) return -1;
1367 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1368 if (!str) return -1;
1369 return (Int_t)str->GetString().Atoi();
1371 //_____________________________________________________________________________________
1372 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
1374 // Author: marian.ivanov@cern.ch
1376 // Create outlier map for CE study
1378 // Return value - outlyer map
1379 // noutlyersCE - number of outlyers
1380 // minSignal - minimal total Q signal
1381 // cutRMSMin - minimal width of the signal in respect to the median
1382 // cutRMSMax - maximal width of the signal in respect to the median
1383 // cutMaxDistT - maximal deviation from time median per chamber
1385 // Outlyers criteria:
1386 // 0. Exclude masked pads
1387 // 1. Exclude first two rows in IROC and last two rows in OROC
1388 // 2. Exclude edge pads
1389 // 3. Exclude channels with too large variations
1390 // 4. Exclude pads with too small signal
1391 // 5. Exclude signal with outlyers RMS
1392 // 6. Exclude channels to far from the chamber median
1394 //create outlier map
1395 AliTPCCalPad *out=ceOut;
1396 if (!out) out= new AliTPCCalPad("outCE","outCE");
1397 AliTPCCalROC *rocMasked=0x0;
1398 if (!fCETmean) return 0;
1399 if (!fCETrms) return 0;
1400 if (!fCEQmean) return 0;
1402 //loop over all channels
1404 Double_t rmsMedian = fCETrms->GetMedian();
1405 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1406 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1407 if (!rocData) continue;
1408 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1409 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1410 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1411 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1412 Double_t trocMedian = rocData->GetMedian();
1414 if (!rocData || !rocCEQ || !rocCETrms || !rocData) {
1415 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1421 UInt_t nrows=rocData->GetNrows();
1422 for (UInt_t irow=0;irow<nrows;++irow){
1423 UInt_t npads=rocData->GetNPads(irow);
1424 for (UInt_t ipad=0;ipad<npads;++ipad){
1425 rocOut->SetValue(irow,ipad,0);
1426 Float_t valTmean=rocData->GetValue(irow,ipad);
1427 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1428 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1429 //0. exclude masked pads
1430 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1431 rocOut->SetValue(irow,ipad,1);
1434 //1. exclude first two rows in IROC and last two rows in OROC
1436 if (irow<2) rocOut->SetValue(irow,ipad,1);
1438 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1440 //2. exclude edge pads
1441 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1442 //exclude values that are exactly 0
1443 if ( TMath::Abs(valTmean)<kAlmost0) {
1444 rocOut->SetValue(irow,ipad,1);
1447 //3. exclude channels with too large variations
1448 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1449 rocOut->SetValue(irow,ipad,1);
1453 //4. exclude channels with too small signal
1454 if (valQmean<minSignal) {
1455 rocOut->SetValue(irow,ipad,1);
1459 //5. exclude channels with too small rms
1460 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1461 rocOut->SetValue(irow,ipad,1);
1465 //6. exclude channels to far from the chamber median
1466 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1467 rocOut->SetValue(irow,ipad,1);
1478 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1480 // Author: marian.ivanov@cern.ch
1482 // Create outlier map for Pulser
1484 // Return value - outlyer map
1485 // noutlyersPulser - number of outlyers
1486 // cutTime - absolute cut - distance to the median of chamber
1487 // cutnRMSQ - nsigma cut from median q distribution per chamber
1488 // cutnRMSrms - nsigma cut from median rms distribution
1489 // Outlyers criteria:
1490 // 0. Exclude masked pads
1491 // 1. Exclude time outlyers (default 3 time bins)
1492 // 2. Exclude q outlyers (default 5 sigma)
1493 // 3. Exclude rms outlyers (default 5 sigma)
1495 AliTPCCalPad *out=pulserOut;
1496 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1497 AliTPCCalROC *rocMasked=0x0;
1498 if (!fPulserTmean) return 0;
1499 if (!fPulserTrms) return 0;
1500 if (!fPulserQmean) return 0;
1502 //loop over all channels
1504 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1505 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1506 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1507 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1508 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1509 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1511 Double_t rocMedianT = rocData->GetMedian();
1512 Double_t rocMedianQ = rocPulserQ->GetMedian();
1513 Double_t rocRMSQ = rocPulserQ->GetRMS();
1514 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1515 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1516 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1517 rocOut->SetValue(ichannel,0);
1518 Float_t valTmean=rocData->GetValue(ichannel);
1519 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1520 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1521 Float_t valMasked =0;
1522 if (rocMasked) valMasked = rocMasked->GetValue(ichannel);
1524 if (valMasked>0.5) isOut=1;
1525 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1526 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1527 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1528 rocOut->SetValue(ichannel,isOut);
1529 if (isOut) noutliersPulser++;
1536 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1538 // Author : Marian Ivanov
1539 // Create pad time0 correction map using information from the CE and from pulser
1542 // Return PadTime0 to be used for time0 relative alignment
1543 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1544 // using $ALICE_ROOT/TPC/script/gui application
1546 // fitResultsA - fitParameters A side
1547 // fitResultsC - fitParameters C side
1548 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1549 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1553 // 1. Find outlier map for CE
1554 // 2. Find outlier map for Pulser
1555 // 3. Replace outlier by median at given sector (median without outliers)
1556 // 4. Substract from the CE data pulser
1557 // 5. Fit the CE with formula
1558 // 5.1) (IROC-OROC) offset
1562 // 5.5) (IROC-OROC)*(lx-xmid)
1564 // 6. Substract gy fit dependence from the CE data
1565 // 7. Add pulser back to CE data
1566 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1567 // 9. return CE data
1569 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1570 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1573 const char *formulaIn="(-1.+2.*(sector<36))*0.5++gx++gy++(lx-134.)++(-1.+2.*(sector<36))*0.5*(lx-134)++((ly/lx)^2/(0.1763)^2)";
1574 // output for fit formula
1575 const char *formulaAll="1++(-1.+2.*(sector<36))*0.5++gx++gy++(lx-134.)++(-1.+2.*(sector<36))*0.5*(lx-134)++((ly/lx)^2/(0.1763)^2)";
1576 // gy part of formula
1577 const char *formulaOut="0++0*(-1.+2.*(sector<36))*0.5++0*gx++gy++0*(lx-134.)++0*(-1.+2.*(sector<36))*0.5*(lx-134)++0*((ly/lx)^2/(0.1763)^2)";
1580 if (!fCETmean) return 0;
1581 Double_t pgya,pgyc,pchi2a,pchi2c;
1582 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1583 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1585 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1586 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1587 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1588 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1589 padPulser->SetName("padPulser");
1590 padPulserOut->SetName("padPulserOut");
1591 padCE->SetName("padCE");
1592 padCEIn->SetName("padCEIn");
1593 padCEOut->SetName("padCEOut");
1594 padOut->SetName("padOut");
1597 // make combined outlyers map
1598 // and replace outlyers in maps with median for chamber
1600 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1601 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1602 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1603 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1604 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1605 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1606 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1607 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1608 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1609 if (rocPulserOut->GetValue(ichannel)>0) {
1610 rocPulser->SetValue(ichannel,pulserMedian);
1611 rocOut->SetValue(ichannel,1);
1613 if (rocCEOut->GetValue(ichannel)>0) {
1614 rocCE->SetValue(ichannel,ceMedian);
1615 rocOut->SetValue(ichannel,1);
1620 // remove pulser time 0
1622 padCE->Add(padPulser,-1);
1627 Float_t chi2Af,chi2Cf;
1628 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1632 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1633 padCEFitGY->SetName("padCEFitGy");
1635 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1636 padCEFit->SetName("padCEFit");
1638 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1639 padCEDiff->SetName("padCEDiff");
1640 padCEDiff->Add(padCEFit,-1.);
1643 padCE->Add(padCEFitGY,-1.);
1645 padCE->Add(padPulser,1.);
1646 Double_t padmedian = padCE->GetMedian();
1647 padCE->Add(-padmedian); // normalize to median
1649 // Replace outliers by fit value - median of diff per given chamber -GY fit
1651 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1652 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1653 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1654 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1655 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1656 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1658 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1659 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1660 if (rocOut->GetValue(ichannel)==0) continue;
1661 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1662 rocCE->SetValue(ichannel,value);
1668 //dump to the file - result can be visualized
1669 AliTPCPreprocessorOnline preprocesor;
1670 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1671 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1672 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1673 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1675 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1676 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1678 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1679 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1680 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1681 preprocesor.DumpToFile(dumpfile);
1684 delete padPulserOut;
1697 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1699 // find the closest point to xref in x direction
1700 // return dx and value
1704 if(!graph) return 0;
1705 if(graph->GetN() < 1) return 0;
1708 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1709 if (index<0) index=0;
1710 if(graph->GetN()==1) {
1711 dx = xref-graph->GetX()[index];
1714 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1715 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1716 dx = xref-graph->GetX()[index];
1718 y = graph->GetY()[index];
1722 Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1724 // Get the correction of the trigger offset
1725 // combining information from the laser track calibration
1726 // and from cosmic calibration
1729 // timeStamp - tim stamp in seconds
1730 // deltaT - integration period to calculate offset
1731 // deltaTLaser -max validity of laser data
1732 // valType - 0 - median, 1- mean
1734 // Integration vaues are just recomendation - if not possible to get points
1735 // automatically increase the validity by factor 2
1736 // (recursive algorithm until one month of data taking)
1739 const Float_t kLaserCut=0.0005;
1740 const Int_t kMaxPeriod=3600*24*30*12; // one year max
1741 const Int_t kMinPoints=20;
1743 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1745 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1747 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1748 if (!array) return 0;
1750 TGraphErrors *laserA[3]={0,0,0};
1751 TGraphErrors *laserC[3]={0,0,0};
1752 TGraphErrors *cosmicAll=0;
1753 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1754 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1755 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1758 if (!cosmicAll) return 0;
1759 Int_t nmeasC=cosmicAll->GetN();
1760 Float_t *tdelta = new Float_t[nmeasC];
1762 for (Int_t i=0;i<nmeasC;i++){
1763 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1764 Float_t ccosmic=cosmicAll->GetY()[i];
1765 Double_t yA=0,yC=0,dA=0,dC=0;
1766 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1767 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1768 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1769 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1771 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1773 if (TMath::Abs(yA-yC)<kLaserCut) {
1776 if (i%2==0) claser=yA;
1777 if (i%2==1) claser=yC;
1779 tdelta[nused]=ccosmic-claser;
1782 if (nused<kMinPoints &&deltaT<kMaxPeriod) {
1784 return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1786 if (nused<kMinPoints) {
1788 //AliWarning("AliFatal: No time offset calibration available\n");
1791 Double_t median = TMath::Median(nused,tdelta);
1792 Double_t mean = TMath::Mean(nused,tdelta);
1794 return (valType==0) ? median:mean;
1797 Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1799 // Get the correction of the drift velocity
1800 // combining information from the laser track calibration
1801 // and from cosmic calibration
1803 // dist - return value - distance to closest point in graph
1805 // timeStamp - tim stamp in seconds
1806 // deltaT - integration period to calculate time0 offset
1807 // deltaTLaser -max validity of laser data
1808 // valType - 0 - median, 1- mean
1810 // Integration vaues are just recomendation - if not possible to get points
1811 // automatically increase the validity by factor 2
1812 // (recursive algorithm until one month of data taking)
1816 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1818 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1820 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1821 if (!array) return 0;
1822 TGraphErrors *cosmicAll=0;
1823 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1824 if (!cosmicAll) return 0;
1826 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1828 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1829 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1830 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1831 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1838 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1839 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1840 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1842 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1843 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1844 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1845 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1847 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1848 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1854 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1857 // Create a default name for the gui file
1860 return Form("guiRefTreeRun%s.root",GetRefValidity());
1863 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1866 // Create a gui reference tree
1867 // if dirname and filename are empty default values will be used
1868 // this is the recommended way of using this function
1869 // it allows to check whether a file with the given run validity alredy exists
1871 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1872 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1876 TString file=filename;
1877 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1879 AliTPCPreprocessorOnline prep;
1880 //noise and pedestals
1881 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1882 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1883 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1885 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1886 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1887 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1888 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1890 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1891 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1892 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1893 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1895 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1896 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1897 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1898 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1899 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1901 AliTPCdataQA *dataQA=fRefDataQA;
1903 if (dataQA->GetNLocalMaxima())
1904 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1905 if (dataQA->GetMaxCharge())
1906 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1907 if (dataQA->GetMeanCharge())
1908 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1909 if (dataQA->GetNoThreshold())
1910 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1911 if (dataQA->GetNTimeBins())
1912 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1913 if (dataQA->GetNPads())
1914 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1915 if (dataQA->GetTimePosition())
1916 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1918 prep.DumpToFile(file.Data());
1922 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1924 // Get the correction of the drift velocity using the offline laser tracks calbration
1927 // timeStamp - tim stamp in seconds
1928 // deltaT - integration period to calculate time0 offset
1929 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1930 // Note in case no data form both A and C side - the value from active side used
1931 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1933 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1936 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1938 // Get the correction of the drift velocity using the online laser tracks calbration
1941 // timeStamp - tim stamp in seconds
1942 // deltaT - integration period to calculate time0 offset
1943 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1944 // Note in case no data form both A and C side - the value from active side used
1945 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1947 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1948 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1949 if (!param) return 0;
1951 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1952 dv*=param->GetDriftV()/2.61301900000000000e+06;
1953 if (dv>1e-20) dv=1/dv-1;
1956 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1958 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
1959 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
1965 AliTPCCalibVdrift corr(temp,press,0);
1966 corrPTA=corr.GetPTRelative(timeStamp,0);
1967 corrPTC=corr.GetPTRelative(timeStamp,1);
1970 if (side==0) dv -= corrPTA;
1971 if (side==1) dv -= corrPTC;
1972 if (side==2) dv -= (corrPTA+corrPTC)/2;
1977 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1978 Int_t side, TObjArray * const array){
1980 // common drift velocity retrieval for online and offline method
1982 TGraphErrors *grlaserA=0;
1983 TGraphErrors *grlaserC=0;
1984 Double_t vlaserA=0, vlaserC=0;
1985 if (!array) return 0;
1986 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1987 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1989 if (grlaserA && grlaserA->GetN()>0) {
1990 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1991 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1992 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1994 if (grlaserC && grlaserC->GetN()>0) {
1995 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1996 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1997 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1999 if (side==0) return vlaserA;
2000 if (side==1) return vlaserC;
2001 Double_t mdrift=(vlaserA+vlaserC)*0.5;
2002 if (!grlaserA) return vlaserC;
2003 if (!grlaserC) return vlaserA;
2008 Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
2010 // Get the correction of the drift velocity using the CE laser data
2011 // combining information from the CE, laser track calibration
2012 // and P/T calibration
2015 // timeStamp - tim stamp in seconds
2016 // deltaT - integration period to calculate time0 offset
2017 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2018 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
2019 if (!arrT) return 0;
2020 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2021 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2022 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2025 Double_t corrPTA = 0, corrPTC=0;
2026 Double_t ltime0A = 0, ltime0C=0;
2028 Double_t corrA=0, corrC=0;
2029 Double_t timeA=0, timeC=0;
2030 const Double_t kEpsilon = 0.00001;
2031 TGraph *graphA = (TGraph*)arrT->At(72);
2032 TGraph *graphC = (TGraph*)arrT->At(73);
2033 if (!graphA && !graphC) return 0.;
2034 if (graphA &&graphA->GetN()>0) {
2035 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2036 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2037 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2038 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2039 if(ltime0A < kEpsilon) return 0;
2040 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
2041 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2044 if (graphC&&graphC->GetN()>0){
2045 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2046 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2047 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2048 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2049 if(ltime0C < kEpsilon) return 0;
2050 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
2051 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2055 if (side ==0 ) return corrA;
2056 if (side ==1 ) return corrC;
2057 Double_t corrM= (corrA+corrC)*0.5;
2058 if (!graphA) corrM=corrC;
2059 if (!graphC) corrM=corrA;
2063 Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2065 // return drift velocity using the TPC-ITS matchin method
2066 // return also distance to the closest point
2068 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2069 TGraphErrors *graph=0;
2071 if (!array) return 0;
2073 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2074 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_TOFB_TPC_DRIFTVD");
2075 if (!graph) return 0;
2077 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2078 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2082 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2084 // Get time dependent time 0 (trigger delay in cm) correction
2086 // timestamp - timestamp
2089 // Notice - Extrapolation outside of calibration range - using constant function
2091 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2092 TGraphErrors *graph=0;
2094 if (!array) return 0;
2095 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2096 if (!graph) graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_T0");
2097 if (!graph) return 0;
2099 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2100 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2108 Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2110 // VERY obscure method - we need something in framework
2111 // Find the TPC runs with temperature OCDB entry
2112 // cache the start and end of the run
2114 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2115 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2116 if (!storage) return 0;
2117 TString path=storage->GetURI();
2121 if (path.Contains("local")){ // find the list if local system
2122 path.ReplaceAll("local://","");
2123 path+="TPC/Calib/Temperature";
2124 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2126 runsT=gSystem->GetFromPipe(command);
2128 TObjArray *arr= runsT.Tokenize("r");
2131 TArrayI indexes(arr->GetEntries());
2132 TArrayI runs(arr->GetEntries());
2134 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2135 Int_t irunN = atoi(arr->At(irun)->GetName());
2136 if (irunN<startRun) continue;
2137 if (irunN>stopRun) continue;
2138 runs[naccept]=irunN;
2143 fRunsStart.Set(fRuns.fN);
2144 fRunsStop.Set(fRuns.fN);
2145 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2146 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2149 AliCDBEntry * entry = 0;
2150 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2151 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2152 if (!entry) continue;
2153 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2154 if (!tmpRun) continue;
2155 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2156 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2157 //AliInfo(Form("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec()));
2163 Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2165 // binary search - find the run for given time stamp
2167 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2168 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2170 for (Int_t index=index0; index<=index1; index++){
2171 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2173 AliInfo(Form("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime));
2176 if (cindex<0) cindex =(index0+index1)/2;
2180 return fRuns[cindex];
2187 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2189 // filter outlyer measurement
2190 // Only points around median +- sigmaCut filtered
2191 if (!graph) return 0;
2193 Int_t npoints0 = graph->GetN();
2198 if (npoints0<kMinPoints) return 0;
2200 Double_t *outx=new Double_t[npoints0];
2201 Double_t *outy=new Double_t[npoints0];
2202 for (Int_t iter=0; iter<3; iter++){
2204 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2205 if (graph->GetY()[ipoint]==0) continue;
2206 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2207 outx[npoints] = graph->GetX()[ipoint];
2208 outy[npoints] = graph->GetY()[ipoint];
2211 if (npoints<=1) break;
2212 medianY =TMath::Median(npoints,outy);
2213 rmsY =TMath::RMS(npoints,outy);
2216 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2223 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2225 // filter outlyer measurement
2226 // Only points around median +- cut filtered
2227 if (!graph) return 0;
2229 Int_t npoints0 = graph->GetN();
2234 if (npoints0<kMinPoints) return 0;
2236 Double_t *outx=new Double_t[npoints0];
2237 Double_t *outy=new Double_t[npoints0];
2238 for (Int_t iter=0; iter<3; iter++){
2240 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2241 if (graph->GetY()[ipoint]==0) continue;
2242 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2243 outx[npoints] = graph->GetX()[ipoint];
2244 outy[npoints] = graph->GetY()[ipoint];
2247 if (npoints<=1) break;
2248 medianY =TMath::Median(npoints,outy);
2249 // rmsY =TMath::RMS(npoints,outy);
2252 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2260 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2262 // filter outlyer measurement
2263 // Only points with normalized errors median +- sigmaCut filtered
2265 Int_t kMinPoints=10;
2266 Int_t npoints0 = graph->GetN();
2268 Float_t medianErr=0, rmsErr=0;
2271 if (npoints0<kMinPoints) return 0;
2273 Double_t *outx=new Double_t[npoints0];
2274 Double_t *outy=new Double_t[npoints0];
2275 Double_t *erry=new Double_t[npoints0];
2276 Double_t *nerry=new Double_t[npoints0];
2277 Double_t *errx=new Double_t[npoints0];
2279 for (Int_t iter=0; iter<3; iter++){
2281 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2282 nerry[npoints] = graph->GetErrorY(ipoint);
2283 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2284 erry[npoints] = graph->GetErrorY(ipoint);
2285 outx[npoints] = graph->GetX()[ipoint];
2286 outy[npoints] = graph->GetY()[ipoint];
2287 errx[npoints] = graph->GetErrorY(ipoint);
2290 if (npoints==0) break;
2291 medianErr=TMath::Median(npoints,erry);
2292 medianY =TMath::Median(npoints,outy);
2293 rmsErr =TMath::RMS(npoints,erry);
2295 TGraphErrors *graphOut=0;
2296 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2306 void AliTPCcalibDButil::Sort(TGraph *graph){
2308 // sort array - neccessay for approx
2310 Int_t npoints = graph->GetN();
2311 Int_t *indexes=new Int_t[npoints];
2312 Double_t *outx=new Double_t[npoints];
2313 Double_t *outy=new Double_t[npoints];
2314 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2315 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2316 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2317 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2318 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2324 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2326 // smmoth graph - mean on the interval
2329 Int_t npoints = graph->GetN();
2330 Double_t *outy=new Double_t[npoints];
2332 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2333 Double_t lx=graph->GetX()[ipoint];
2334 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2335 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2336 if (index0<0) index0=0;
2337 if (index1>=npoints-1) index1=npoints-1;
2338 if ((index1-index0)>1){
2339 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2341 outy[ipoint]=graph->GetY()[ipoint];
2344 // TLinearFitter fitter(3,"pol2");
2345 // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2346 // Double_t lx=graph->GetX()[ipoint];
2347 // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2348 // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2349 // if (index0<0) index0=0;
2350 // if (index1>=npoints-1) index1=npoints-1;
2351 // fitter.ClearPoints();
2352 // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2353 // if ((index1-index0)>1){
2354 // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2356 // outy[ipoint]=graph->GetY()[ipoint];
2362 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2363 graph->GetY()[ipoint] = outy[ipoint];
2368 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2370 // Use constant interpolation outside of range
2373 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2377 if (graph->GetN()<1){
2378 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
2383 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2384 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2386 // AliInfo(Form("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref)));
2388 if(graph->GetN()==1)
2389 return graph->Eval(graph->GetX()[0]);
2392 return graph->Eval(xref);
2395 Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2397 // Use constant interpolation outside of range also for spline fits
2400 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2403 if (graph->GetKnots()<1){
2404 AliInfoGeneral("AliTPCcalibDButil","AliTPCcalibDButil::EvalGraphConst: Empty graph");
2407 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2408 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2409 return graph->Eval( xref);
2412 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2414 // Filter DCS sensor information
2415 // ymin - minimal value
2417 // maxdy - maximal deirivative
2418 // sigmaCut - cut on values and derivative in terms of RMS distribution
2419 // Return value - accepted fraction
2423 // 0. Calculate median and rms of values in specified range
2424 // 1. Filter out outliers - median+-sigmaCut*rms
2425 // values replaced by median
2427 AliSplineFit * fit = sensor->GetFit();
2428 if (!fit) return 0.;
2429 Int_t nknots = fit->GetKnots();
2436 Double_t *yin0 = new Double_t[nknots];
2437 Double_t *yin1 = new Double_t[nknots];
2440 for (Int_t iknot=0; iknot< nknots; iknot++){
2441 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2442 yin0[naccept] = fit->GetY0()[iknot];
2443 yin1[naccept] = fit->GetY1()[iknot];
2444 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2456 Double_t medianY0=0, medianY1=0;
2457 Double_t rmsY0 =0, rmsY1=0;
2458 medianY0 = TMath::Median(naccept, yin0);
2459 medianY1 = TMath::Median(naccept, yin1);
2460 rmsY0 = TMath::RMS(naccept, yin0);
2461 rmsY1 = TMath::RMS(naccept, yin1);
2464 // 1. Filter out outliers - median+-sigmaCut*rms
2465 // values replaced by median
2466 // if replaced the derivative set to 0
2468 for (Int_t iknot=0; iknot< nknots; iknot++){
2470 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2471 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2472 if (nknots<2) fit->GetY1()[iknot]=0;
2473 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2475 fit->GetY0()[iknot]=medianY0;
2476 fit->GetY1()[iknot]=0;
2483 return Float_t(naccept)/Float_t(nknots);
2486 Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2488 // Filter temperature array
2489 // tempArray - array of temperatures -
2490 // ymin - minimal accepted temperature - default 15
2491 // ymax - maximal accepted temperature - default 22
2492 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2493 // return value - fraction of filtered sensors
2494 const Double_t kMaxDy=0.1;
2495 Int_t nsensors=tempArray->NumSensors();
2496 if (nsensors==0) return 0.;
2498 for (Int_t isensor=0; isensor<nsensors; isensor++){
2499 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2500 if (!sensor) continue;
2501 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2502 if (sensor->GetFit()==0){
2504 tempArray->RemoveSensorNum(isensor);
2509 return Float_t(naccept)/Float_t(nsensors);
2513 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2516 // Input parameters:
2517 // deltaT - smoothing window (in seconds)
2518 // cutAbs - max distance of the time info to the median (in time bins)
2519 // cutSigma - max distance (in the RMS)
2520 // pcstream - optional debug streamer to store original and filtered info
2521 // Hardwired parameters:
2522 // kMinPoints =10; // minimal number of points to define the CE
2523 // kMinSectors=12; // minimal number of sectors to define sideCE
2525 // 0. Filter almost emty graphs (kMinPoints=10)
2526 // 1. calculate median and RMS per side
2527 // 2. Filter graphs - in respect with side medians
2528 // - cutAbs and cutDelta used
2529 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2530 // 4. Calculate mean for A side and C side
2532 const Int_t kMinPoints =10; // minimal number of points to define the CE
2533 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2534 const Int_t kMinTime =400; // minimal arrival time of CE
2535 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2537 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2538 if (!cearray) return;
2543 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2544 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2545 if ( tempMapCE && cavernPressureCE){
2547 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2548 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2549 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2552 // recalculate P/T correction map for time of the CE
2553 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2554 driftCalib->SetName("driftPTCE");
2555 driftCalib->SetTitle("driftPTCE");
2556 cearray->AddLast(driftCalib);
2560 // 0. Filter almost emty graphs
2563 for (Int_t i=0; i<72;i++){
2564 TGraph *graph= (TGraph*)arrT->At(i);
2565 if (!graph) continue;
2567 if (graph->GetN()<kMinPoints){
2569 delete graph; // delete empty graph
2572 if (tmin<0) tmin = graph->GetX()[0];
2573 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2575 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2576 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2579 // 1. calculate median and RMS per side
2581 TArrayF arrA(100000), arrC(100000);
2583 Double_t medianA=0, medianC=0;
2584 Double_t rmsA=0, rmsC=0;
2585 for (Int_t isec=0; isec<72;isec++){
2586 TGraph *graph= (TGraph*)arrT->At(isec);
2587 if (!graph) continue;
2588 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2589 if (graph->GetY()[ipoint]<kMinTime) continue;
2590 if (nA>=arrA.fN) arrA.Set(nA*2);
2591 if (nC>=arrC.fN) arrC.Set(nC*2);
2592 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2593 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2597 medianA=TMath::Median(nA,arrA.fArray);
2598 rmsA =TMath::RMS(nA,arrA.fArray);
2601 medianC=TMath::Median(nC,arrC.fArray);
2602 rmsC =TMath::RMS(nC,arrC.fArray);
2605 // 2. Filter graphs - in respect with side medians
2607 TArrayD vecX(100000), vecY(100000);
2608 for (Int_t isec=0; isec<72;isec++){
2609 TGraph *graph= (TGraph*)arrT->At(isec);
2610 if (!graph) continue;
2611 Double_t median = (isec%36<18) ? medianA: medianC;
2612 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2614 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2615 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
2616 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2617 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2618 vecX[naccept]= graph->GetX()[ipoint];
2619 vecY[naccept]= graph->GetY()[ipoint];
2622 if (naccept<kMinPoints){
2623 arrT->AddAt(0,isec);
2624 delete graph; // delete empty graph
2627 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2629 arrT->AddAt(graph2,isec);
2632 // 3. Cut in respect wit the graph median
2634 for (Int_t i=0; i<72;i++){
2635 TGraph *graph= (TGraph*)arrT->At(i);
2636 if (!graph) continue;
2640 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2641 if (!graphTS0) continue;
2642 if (graphTS0->GetN()<kMinPoints) {
2648 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2649 if (!graphTS) continue;
2651 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2653 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2654 (*pcstream)<<"filterCE"<<
2659 "graphTS0.="<<graphTS0<<
2660 "graphTS.="<<graphTS<<
2664 arrT->AddAt(graphTS,i);
2668 // Recalculate the mean time A side C side
2670 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2671 Int_t meanPoints=(nA+nC)/72; // mean number of points
2672 for (Int_t itime=0; itime<200; itime++){
2674 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2675 for (Int_t i=0; i<72;i++){
2676 TGraph *graph= (TGraph*)arrT->At(i);
2677 if (!graph) continue;
2678 if (graph->GetN()<(meanPoints/4)) continue;
2679 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2680 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2684 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2685 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2686 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2687 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2690 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2691 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2693 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2694 (*pcstream)<<"filterAC"<<
2703 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2704 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2705 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2706 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2707 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2708 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2711 if (nA<kMinSectors) arrT->AddAt(0,72);
2712 else arrT->AddAt(graphTSA,72);
2713 if (nC<kMinSectors) arrT->AddAt(0,73);
2714 else arrT->AddAt(graphTSC,73);
2718 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2720 // Filter Drift velocity measurement using the tracks
2721 // 0. remove outlyers - error based
2725 const Int_t kMinPoints=1; // minimal number of points to define value
2726 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2729 for (Int_t i=0; i<arrT->GetEntries();i++){
2730 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
2731 if (!graph) continue;
2732 if (graph->GetN()<kMinPoints){
2737 TGraphErrors *graph2 = NULL;
2738 if(graph->GetN()<10) {
2739 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2741 delete graph; arrT->AddAt(0,i); continue;
2745 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2747 delete graph; arrT->AddAt(0,i); continue;
2750 if (graph2->GetN()<1) {
2751 delete graph; arrT->AddAt(0,i); continue;
2753 graph2->SetName(graph->GetName());
2754 graph2->SetTitle(graph->GetTitle());
2755 arrT->AddAt(graph2,i);
2757 (*pcstream)<<"filterTracks"<<
2762 "graph2.="<<graph2<<
2773 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2776 // get laser time offset
2777 // median around timeStamp+-deltaT
2778 // QA - chi2 needed for later usage - to be added
2779 // - currently cut on error
2782 Double_t kMinDelay=0.01;
2783 Double_t kMinDelayErr=0.0001;
2785 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2786 if (!array) return 0;
2787 TGraphErrors *tlaser=0;
2789 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2790 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2792 if (!tlaser) return 0;
2793 Int_t npoints0= tlaser->GetN();
2794 if (npoints0==0) return 0;
2795 Double_t *xlaser = new Double_t[npoints0];
2796 Double_t *ylaser = new Double_t[npoints0];
2798 for (Int_t i=0;i<npoints0;i++){
2800 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2801 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2802 xlaser[npoints]=tlaser->GetX()[npoints];
2803 ylaser[npoints]=tlaser->GetY()[npoints];
2808 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2809 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2810 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2811 if (index0<0) index0=0;
2812 if (index1>=npoints-1) index1=npoints-1;
2813 if (index1-index0<kMinPoints) {
2819 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2820 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2829 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2831 // Filter Goofie data
2832 // goofieArray - points will be filtered
2833 // deltaT - smmothing time window
2834 // cutSigma - outler sigma cut in rms
2835 // minVn, maxVd- range absolute cut for variable vd/pt
2838 // Ignore goofie if not enough points
2840 const Int_t kMinPoints = 3;
2843 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2844 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2845 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2846 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2847 if (!graphvd) return;
2848 if (graphvd->GetN()<kMinPoints){
2850 goofieArray->GetSensorNum(2)->SetGraph(0);
2854 // 1. Caluclate medians of critical variables
2860 Double_t medianpt=0;
2861 Double_t medianvd=0, sigmavd=0;
2862 Double_t medianan=0;
2863 Double_t medianaf=0;
2864 Int_t entries=graphvd->GetN();
2865 Double_t yvdn[10000];
2868 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2869 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2870 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2871 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2872 yvdn[nvd++]=graphvd->GetY()[ipoint];
2874 if (nvd<kMinPoints){
2876 goofieArray->GetSensorNum(2)->SetGraph(0);
2880 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2881 if (nuni>=kMinPoints){
2882 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2884 medianvd = TMath::Median(nvd, yvdn);
2887 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2888 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2889 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2890 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2891 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2892 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2900 // 2. Make outlyer graph
2903 TGraph graphOut(*graphvd);
2904 for (Int_t i=0; i<entries;i++){
2906 Bool_t isOut=kFALSE;
2907 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2908 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2910 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2912 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2913 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2914 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2915 graphOut.GetY()[i]= (isOut)?1:0;
2918 if (nOK<kMinPoints) {
2920 goofieArray->GetSensorNum(2)->SetGraph(0);
2924 // 3. Filter out outlyers - and smooth
2926 TVectorF vmedianArray(goofieArray->NumSensors());
2927 TVectorF vrmsArray(goofieArray->NumSensors());
2928 Double_t xnew[10000];
2929 Double_t ynew[10000];
2931 junk.SetOwner(kTRUE);
2935 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2937 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2938 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2940 if (!sensor) continue;
2941 graphOld = sensor->GetGraph();
2943 sensor->SetGraph(0);
2945 for (Int_t i=0;i<entries;i++){
2946 if (graphOut.GetY()[i]>0.5) continue;
2947 xnew[nused]=graphOld->GetX()[i];
2948 ynew[nused]=graphOld->GetY()[i];
2951 graphNew = new TGraph(nused,xnew,ynew);
2952 junk.AddLast(graphNew);
2953 junk.AddLast(graphOld);
2955 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2957 junk.AddLast(graphNew0);
2958 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2960 junk.AddLast(graphNew1);
2961 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2963 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2964 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2965 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2966 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2967 // AliInfo(Form("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]));
2968 vmedianArray[isensor]=median;
2974 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2975 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2976 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2977 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2978 (*pcstream)<<"goofieA"<<
2979 Form("isOK_%d.=",isensor)<<isOK<<
2980 Form("s_%d.=",isensor)<<sensor<<
2981 Form("gr_%d.=",isensor)<<graphOld<<
2982 Form("gr0_%d.=",isensor)<<graphNew0<<
2983 Form("gr1_%d.=",isensor)<<graphNew1<<
2984 Form("gr2_%d.=",isensor)<<graphNew2;
2985 if (isOK) sensor->SetGraph(graphNew2);
2987 (*pcstream)<<"goofieA"<<
2988 "vmed.="<<&vmedianArray<<
2989 "vrms.="<<&vrmsArray<<
2991 junk.Delete(); // delete temoprary graphs
2999 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
3001 // Make a statistic matrix
3002 // Input parameters:
3003 // array - TObjArray of AliRelKalmanAlign
3004 // minFraction - minimal ration of accepted tracks
3005 // minStat - minimal statistic (number of accepted tracks)
3006 // maxvd - maximal deviation for the 1
3008 // columns - Mean, Median, RMS
3009 // row - parameter type (rotation[3], translation[3], drift[3])
3010 if (!array) return 0;
3011 if (array->GetEntries()<=0) return 0;
3012 // Int_t entries = array->GetEntries();
3013 Int_t entriesFast = array->GetEntriesFast();
3015 TVectorD *valArray[9];
3016 for (Int_t i=0; i<9; i++){
3017 valArray[i] = new TVectorD(entriesFast);
3020 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3021 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3022 if (!kalman) continue;
3023 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3024 if (kalman->GetNUpdates()<minStat) continue;
3025 if (Float_t(kalman->GetNUpdates())/Float_t(kalman->GetNTracks())<minFraction) continue;
3026 kalman->GetState(state);
3027 for (Int_t ipar=0; ipar<9; ipar++)
3028 (*valArray[ipar])[naccept]=state[ipar];
3031 //if (naccept<2) return 0;
3032 if (naccept<1) return 0;
3033 TMatrixD *pstat=new TMatrixD(9,3);
3034 TMatrixD &stat=*pstat;
3035 for (Int_t ipar=0; ipar<9; ipar++){
3036 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3037 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3038 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3040 for (Int_t i=0; i<9; i++){
3047 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
3049 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3051 // array - input array
3052 // stat - mean parameters statistic
3054 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3055 if (!array) return 0;
3056 if (array->GetEntries()<=0) return 0;
3057 if (!(&stat)) return 0;
3058 // error increase in 1 hour
3059 const Double_t kerrsTime[9]={
3060 0.00001, 0.00001, 0.00001,
3061 0.001, 0.001, 0.001,
3062 0.002, 0.01, 0.001};
3065 Int_t entries = array->GetEntriesFast();
3066 TObjArray *sArray= new TObjArray(entries);
3067 AliRelAlignerKalman * sKalman =0;
3069 for (Int_t i=0; i<entries; i++){
3070 Int_t index=(direction)? entries-i-1:i;
3071 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3072 if (!kalman) continue;
3074 kalman->GetState(state);
3075 for (Int_t ipar=0; ipar<9; ipar++){
3076 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3078 if (!sKalman &&isOK) {
3079 sKalman=new AliRelAlignerKalman(*kalman);
3080 sKalman->SetRejectOutliers(kFALSE);
3081 sKalman->SetRunNumber(kalman->GetRunNumber());
3082 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3084 if (!sKalman) continue;
3085 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
3086 for (Int_t ipar=0; ipar<9; ipar++){
3087 // (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3088 // (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3089 // (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3090 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3092 sKalman->SetRunNumber(kalman->GetRunNumber());
3093 if (!isOK) sKalman->SetRunNumber(0);
3094 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3095 if (!isOK) continue;
3096 sKalman->SetRejectOutliers(kFALSE);
3097 sKalman->SetRunNumber(kalman->GetRunNumber());
3098 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3099 sKalman->Merge(kalman);
3100 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3106 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
3108 // Merge 2 RelKalman arrays
3110 // arrayP - rel kalman in direction plus
3111 // arrayM - rel kalman in direction minus
3112 if (!arrayP) return 0;
3113 if (arrayP->GetEntries()<=0) return 0;
3114 if (!arrayM) return 0;
3115 if (arrayM->GetEntries()<=0) return 0;
3117 Int_t entries = arrayP->GetEntriesFast();
3118 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3120 for (Int_t i=0; i<entries; i++){
3121 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3122 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3123 if (!kalmanP) continue;
3124 if (!kalmanM) continue;
3126 AliRelAlignerKalman *kalman = NULL;
3127 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3128 kalman = new AliRelAlignerKalman(*kalmanP);
3129 kalman->Merge(kalmanM);
3131 else if (kalmanP->GetRunNumber() == 0) {
3132 kalman = new AliRelAlignerKalman(*kalmanM);
3134 else if (kalmanM->GetRunNumber() == 0) {
3135 kalman = new AliRelAlignerKalman(*kalmanP);
3140 array->AddAt(kalman,i);
3146 //_____________________________________________________________________________________
3147 TTree* AliTPCcalibDButil::ConnectGainTrees(TString baseDir)
3150 // baseDir: Base directory with the raw Kr calibration trees
3151 // and the trees from the calibQA
3152 // it assumes to following structure below:
3153 // KryptonCalib/<year>/calibKr/calibKr.<year>.<id>.root
3154 // calibQAdEdx/<year>/calibQA.<year>.<perid>.tree.root
3155 // map/treeMapping.root
3159 // === add main tree, which will be a mapping file ================
3160 TFile *fin = TFile::Open(Form("%s/map/treeMapping.root",baseDir.Data()));
3162 TTree *tMain = (TTree*)fin->Get("calPads");
3163 tMain->SetAlias("Pad0","MapPadLength.fElements==7.5"); // pad types
3164 tMain->SetAlias("Pad1","MapPadLength.fElements==10.0");
3165 tMain->SetAlias("Pad2","MapPadLength.fElements==15.0");
3167 tMain->SetAlias("dRowNorm0","(row.fElements/32-1)"); // normalized distance to the center of the chamber in lx (-1,1)
3168 tMain->SetAlias("dRowNorm1","(row.fElements/32-1)"); //
3169 tMain->SetAlias("dRowNorm2","((row.fElements-63)/16-1)"); //
3170 tMain->SetAlias("dRowNorm","(row.fElements<63)*(row.fElements/32-1)+(row.fElements>63)*((row.fElements-63)/16-1)");
3171 tMain->SetAlias("dPhiNorm","(ly.fElements/lx.fElements)/(pi/18)"); // normalized distance to the center of the chamber in phi (-1,1)
3173 tMain->SetAlias("fsector","(sector%36)+(0.5+9*atan2(ly.fElements,lx.fElements)/pi)"); // float sector number
3174 tMain->SetAlias("posEdge","((pi/18.)-abs(atan(ly.fElements/lx.fElements)))*lx.fElements"); //distance to the edge
3176 // === add the krypton calibration trees ==========================
3177 TString inputTreesKrCalib = gSystem->GetFromPipe(Form("ls %s/KryptonCalib/20*/calibKr/*.tree.root",baseDir.Data()));
3178 TObjArray *arrInputTreesKrCalib = inputTreesKrCalib.Tokenize("\n");
3180 for (Int_t itree=0; itree<arrInputTreesKrCalib->GetEntriesFast(); ++itree) {
3181 TFile *fin2 = TFile::Open(arrInputTreesKrCalib->At(itree)->GetName());
3182 TTree *tin = (TTree*)fin2->Get("calPads");
3184 TString friendName=gSystem->BaseName(arrInputTreesKrCalib->At(itree)->GetName());
3185 friendName.ReplaceAll("calibKr.","");
3186 friendName.ReplaceAll(".tree.root","");
3187 friendName="Kr."+friendName;
3188 tMain->AddFriend(tin,friendName.Data());
3192 // TODO: finish implementation of alias via lists
3193 // const Int_t nbranchAlias = 2;
3194 // const char* branchNames[nbranchAlias]={"spectrMean","fitMean"};
3195 // const Int_t nbranchStat = 2;
3196 // const char* statNames[nbranchStat] = {"Median","LTM"};
3198 // for (Int_t iname=0; iname<nbranchAlias; ++iname) {
3199 // TString branchName = TString::Format("%s.%s", friendName.Data(), bNames[i]);
3201 // for (Int_t istat=0; istat<nbranchStat; ++istat) {
3206 tMain->SetAlias((friendName+".spectrMean_LTMRatio").Data(),
3207 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_LTM+0))",
3208 friendName.Data(),friendName.Data()).Data());
3210 tMain->SetAlias((friendName+".spectrMean_MedianRatio").Data(),
3211 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Median+0))",
3212 friendName.Data(),friendName.Data()).Data());
3213 tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
3214 TString::Format("(%s.spectrMean.fElements/(%s.spectrMean_Mean+0))",
3215 friendName.Data(),friendName.Data()).Data());
3216 tMain->SetAlias((friendName+".spectrMean_MeanRatio").Data(),
3217 TString::Format("(%s.spectrMean.fElements/%s.spectrMean_Mean)",
3218 friendName.Data(),friendName.Data()).Data());
3220 tMain->SetAlias((friendName+".fitMean_LTMRatio").Data(),
3221 TString::Format("(%s.fitMean.fElements/(%s.fitMean_LTM+0))",
3222 friendName.Data(),friendName.Data()).Data());
3223 tMain->SetAlias((friendName+".fitRMS_LTMRatio").Data(),
3224 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_LTM+0))",
3225 friendName.Data(),friendName.Data()).Data());
3227 tMain->SetAlias((friendName+".fitMean_MedianRatio").Data(),
3228 TString::Format("(%s.fitMean.fElements/(%s.fitMean_Median+0))",
3229 friendName.Data(),friendName.Data()).Data());
3230 tMain->SetAlias((friendName+".fitRMS_MedianRatio").Data(),
3231 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Median+0))",
3232 friendName.Data(),friendName.Data()).Data());
3233 tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
3234 TString::Format("(%s.fitMean.fElements/(%s.fitMean_Mean+0))",
3235 friendName.Data(),friendName.Data()).Data());
3236 tMain->SetAlias((friendName+".fitRMS_MeanRatio").Data(),
3237 TString::Format("(%s.fitRMS.fElements/(%s.fitRMS_Mean+0))",
3238 friendName.Data(),friendName.Data()).Data());
3239 tMain->SetAlias((friendName+".fitMean_MeanRatio").Data(),
3240 TString::Format("(%s.fitMean.fElements/%s.fitMean_Mean)",
3241 friendName.Data(),friendName.Data()).Data());
3245 // === add the calibQA trees ======================================
3246 TString inputTreesQACalib = gSystem->GetFromPipe(Form("ls %s/calibQAdEdx/20*/*.tree.root",baseDir.Data()));
3247 TObjArray *arrInputTreesQACalib = inputTreesQACalib.Tokenize("\n");
3249 for (Int_t itree=0; itree<arrInputTreesQACalib->GetEntriesFast(); ++itree) {
3250 TFile *fin2 = TFile::Open(arrInputTreesQACalib->At(itree)->GetName());
3251 TTree *tin = (TTree*)fin2->Get("calPads");
3253 TString friendName=gSystem->BaseName(arrInputTreesQACalib->At(itree)->GetName());
3254 friendName.ReplaceAll("calibQA.","");
3255 friendName.ReplaceAll(".tree.root","");
3256 friendName="QA."+friendName;
3257 tMain->AddFriend(tin,friendName.Data());
3260 tMain->SetAlias((friendName+".MaxCharge_LTMRatio").Data(),
3261 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_LTM)",
3262 friendName.Data(),friendName.Data()).Data());
3264 tMain->SetAlias((friendName+".MaxCharge_MedianRatio").Data(),
3265 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Median)",
3266 friendName.Data(),friendName.Data()).Data());
3267 tMain->SetAlias((friendName+".MaxCharge_MeanRatio").Data(),
3268 TString::Format("(%s.MaxCharge.fElements/%s.MaxCharge_Mean)",
3269 friendName.Data(),friendName.Data()).Data());
3271 tMain->SetAlias((friendName+".MeanCharge_LTMRatio").Data(),
3272 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_LTM)",
3273 friendName.Data(),friendName.Data()).Data());
3275 tMain->SetAlias((friendName+".MeanCharge_MedianRatio").Data(),
3276 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Median)",
3277 friendName.Data(),friendName.Data()).Data());
3278 tMain->SetAlias((friendName+".MeanCharge_MeanRatio").Data(),
3279 TString::Format("(%s.MeanCharge.fElements/%s.MeanCharge_Mean)",
3280 friendName.Data(),friendName.Data()).Data());
3288 //_____________________________________________________________________________________
3289 TTree* AliTPCcalibDButil::ConnectPulserTrees(TString baseDir, TTree *tMain)
3292 // baseDir: Base directory with Pulser information
3293 // TTrees are added to the base tree as a friend tree
3295 // === add the calibPulser trees ======================================
3296 TString inputTreesPulserCalib = gSystem->GetFromPipe(Form("ls %s/calibPulser/20*/*.tree.root",baseDir.Data()));
3297 TObjArray *arrInputTreesPulserCalib = inputTreesPulserCalib.Tokenize("\n");
3298 for (Int_t itree=0; itree<arrInputTreesPulserCalib->GetEntriesFast(); ++itree) {
3299 TFile *fin2 = TFile::Open(arrInputTreesPulserCalib->At(itree)->GetName());
3300 TTree *tin = (TTree*)fin2->Get("calPads");
3302 TString friendName=gSystem->BaseName(arrInputTreesPulserCalib->At(itree)->GetName());
3303 friendName.ReplaceAll("calibPulser.","");
3304 friendName.ReplaceAll(".tree.root","");
3305 friendName="Pulser."+friendName;
3306 tMain->AddFriend(tin,friendName.Data());
3309 tMain->SetAlias((friendName+".CEQmean_LTMRatio").Data(),
3310 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_LTM)",
3311 friendName.Data(),friendName.Data()).Data());
3312 tMain->SetAlias((friendName+".CEQmean_MedianRatio").Data(),
3313 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Median)",
3314 friendName.Data(),friendName.Data()).Data());
3315 tMain->SetAlias((friendName+".CEQmean_MeanRatio").Data(),
3316 TString::Format("(%s.CEQmean.fElements/%s.CEQmean_Mean)",
3317 friendName.Data(),friendName.Data()).Data());
3319 tMain->SetAlias((friendName+".CETmean_LTMDelta").Data(),
3320 TString::Format("(%s.CETmean.fElements-%s.CETmean_LTM)",
3321 friendName.Data(),friendName.Data()).Data());
3322 tMain->SetAlias((friendName+".CETmean_MedianDelta").Data(),
3323 TString::Format("(%s.CETmean.fElements-%s.CETmean_Median)",
3324 friendName.Data(),friendName.Data()).Data());
3325 tMain->SetAlias((friendName+".CETmean_MeanDelta").Data(),
3326 TString::Format("(%s.CETmean.fElements-%s.CETmean_Mean)",
3327 friendName.Data(),friendName.Data()).Data());
3333 TTree* AliTPCcalibDButil::ConnectDistortionTrees(TString baseDir, TString selection, TTree *tMain){
3335 // baseDir: Base directory with Distortion information
3336 // TTrees are added to the base tree as a friend tree
3337 // If base tree not provide - first tree from list is used as base
3339 // === add the calibDistortion trees ======================================
3340 //TString inputTreesDistortionCalib = gSystem->GetFromPipe(Form("ls %s/calibDistortion/20*/*.tree.root",baseDir.Data()));
3341 // TString baseDir="$NOTES/reconstruction/distortionFit/"; TTree *tMain=0;
3342 // AliTPCcalibDButil::ConnectDistortionTrees("$NOTES/reconstruction/distortionFit/", "calibTimeResHisto.root", 0);
3344 TString inputTreesDistortionCalib = "";
3345 if (selection.Contains(".list")){
3346 inputTreesDistortionCalib=gSystem->GetFromPipe(Form("cat %s",selection.Data()));
3348 inputTreesDistortionCalib=gSystem->GetFromPipe(Form("find %s -iname \"%s\"",baseDir.Data(),selection.Data()));
3350 TObjArray *arrInputTreesDistortionCalib = inputTreesDistortionCalib.Tokenize("\n");
3352 TString treePrefix="";
3353 for (Int_t itree=0; itree<arrInputTreesDistortionCalib->GetEntriesFast(); ++itree) {
3354 TString fstring = arrInputTreesDistortionCalib->At(itree)->GetName();
3355 if (fstring.Index("Prefix:")==0){ //get tree descriptor
3356 treePrefix=TString(&(fstring[7]));
3359 if (fstring.Index("#")==0){ // ignore comment field
3362 TFile *finput= TFile::Open(arrInputTreesDistortionCalib->At(itree)->GetName());
3363 TString strFile=arrInputTreesDistortionCalib->At(itree)->GetName();
3364 TObjArray *path=strFile.Tokenize("/");
3365 Int_t plength=path->GetEntries();
3366 if (!finput) continue;
3367 TList* list = finput->GetListOfKeys();
3368 Int_t nkeys=list->GetEntries();
3369 for (Int_t ikey=0; ikey<nkeys; ikey++){
3370 TKey * key = (TKey*)list->At(ikey);
3371 if (strstr(key->GetClassName(),"TTree")==0) continue;
3372 TTree * tree = dynamic_cast<TTree*>(finput->Get(list->At(ikey)->GetName()));
3373 if (!tree) continue;
3374 TString friendName="";
3375 if (treePrefix.Length()==0) {
3376 friendName=TString::Format("%s.%s.%s",path->At(plength-3)->GetName(),path->At(plength-2)->GetName(), tree->GetName());
3378 friendName=TString::Format("%s.%s",treePrefix.Data(), tree->GetName());
3380 ::Info("AliTPCcalibDButil::ConnectDistortionTrees","%s",friendName.Data());
3383 tree = dynamic_cast<TTree*>(finput->Get(list->At(ikey)->GetName()));
3385 tMain->AddFriend(tree,friendName.Data());
3388 // tMain->SetAlias("");
3393 //_____________________________________________________________________________________
3394 TTree* AliTPCcalibDButil::ConnectCalPadTrees(TString baseDir, TString pattern, TTree *tMain, Bool_t checkAliases)
3397 // baseDir: Base directory with per Pad information
3398 // TTrees are added to the base tree as a friend tree
3400 // TString baseDir="/hera/alice/fsozzi/summarymaps/calib2/"; // prefix directory with calibration with slash at the end
3401 // TString pattern="QA/*/*root";
3402 // TTree * tree = AliTPCcalibDButil::ConnectCalPadTrees(baseDir,pattern,0); //create tree and attach calibration as friends
3405 // === add the calibPulser trees ======================================
3406 TString inputTreesCalPad = gSystem->GetFromPipe(Form("ls %s/%s",baseDir.Data(), pattern.Data()));
3407 TObjArray *arrInputTreesCalPad = inputTreesCalPad.Tokenize("\n");
3409 for (Int_t itree=0; itree<arrInputTreesCalPad->GetEntriesFast(); ++itree) {
3410 TFile *fin2 = TFile::Open(arrInputTreesCalPad->At(itree)->GetName());
3411 TTree *tin = (TTree*)fin2->Get("calPads");
3413 TString friendName=arrInputTreesCalPad->At(itree)->GetName();
3414 friendName.ReplaceAll("//","/");
3415 friendName.ReplaceAll(baseDir.Data(),"");
3416 friendName.ReplaceAll("^/","");
3417 friendName.ReplaceAll("/",".");
3418 friendName.ReplaceAll(".tree.",".");
3419 friendName.ReplaceAll(".root","");
3420 printf("%s\n", friendName.Data());
3421 ::Info("AliTPCcalibDButil::ConnectCalPadTrees","%s",friendName.Data());
3422 if (tMain==0) tMain=tin;
3423 tMain->AddFriend(tin,friendName.Data());
3424 TObjArray * branches=tin->GetListOfBranches();
3425 Int_t nBranches=branches->GetEntries();
3426 for (Int_t ibranch=0; ibranch<nBranches; ibranch++){
3427 TString bname=branches->At(ibranch)->GetName();
3428 if (bname.Contains(".")>0){
3429 bname.ReplaceAll(".","");
3431 tin->SetAlias((bname).Data(), (bname+".fElements").Data());
3432 tMain->SetAlias((friendName+"."+bname).Data(), (friendName+"."+bname+".fElements").Data());
3434 // make normalized values per chamber
3436 if (branches->FindObject(bname+"_LTM")!=0){
3437 tMain->SetAlias((friendName+"."+bname+"_MeanRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Mean").Data());
3438 tMain->SetAlias((friendName+"."+bname+"_MedianRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_Median").Data());
3439 tMain->SetAlias((friendName+"."+bname+"_LTMRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"_LTM").Data());
3440 tMain->SetAlias((friendName+"."+bname+"_MeanDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Mean").Data());
3441 tMain->SetAlias((friendName+"."+bname+"_MedianDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_Median").Data());
3442 tMain->SetAlias((friendName+"."+bname+"_LTMDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"_LTM").Data());
3444 // if (branches->FindObject(bname+"_Med3.")!=0){
3445 tMain->SetAlias((friendName+"."+bname+"_Med3Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med3.fElements").Data());
3446 tMain->SetAlias((friendName+"."+bname+"_Med5Ratio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Med5.fElements").Data());
3447 tMain->SetAlias((friendName+"."+bname+"_Par6GRatio").Data(), (friendName+"."+bname+".fElements/"+friendName+"."+bname+"Par6G.fElements").Data());
3448 tMain->SetAlias((friendName+"."+bname+"_Med3Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med3.fElements").Data());
3449 tMain->SetAlias((friendName+"."+bname+"_Med5Delta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Med5.fElements").Data());
3450 tMain->SetAlias((friendName+"."+bname+"_Par6GDelta").Data(), (friendName+"."+bname+".fElements-"+friendName+"."+bname+"Par6G.fElements").Data());
3459 // to be implemented