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"
59 const Float_t kAlmost0=1.e-30;
61 ClassImp(AliTPCcalibDButil)
62 AliTPCcalibDButil::AliTPCcalibDButil() :
70 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
82 fRefPedestalMasked(0x0),
86 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
87 fRefPulserMasked(0x0),
94 fRefALTROAcqStart(0x0),
95 fRefALTROAcqStop(0x0),
100 fMapper(new AliTPCmapper(0x0)),
101 fNpulserOutliers(-1),
103 fCETmaxLimitAbs(1.5),
104 fPulTmaxLimitAbs(1.5),
107 fRuns(0), // run list with OCDB info
108 fRunsStart(0), // start time for given run
109 fRunsStop(0) // stop time for given run
115 //_____________________________________________________________________________________
116 AliTPCcalibDButil::~AliTPCcalibDButil()
121 delete fPulserOutlier;
122 delete fRefPulserOutlier;
124 if (fRefPadNoise) delete fRefPadNoise;
125 if (fRefPedestals) delete fRefPedestals;
126 if (fRefPedestalMasked) delete fRefPedestalMasked;
127 if (fRefPulserTmean) delete fRefPulserTmean;
128 if (fRefPulserTrms) delete fRefPulserTrms;
129 if (fRefPulserQmean) delete fRefPulserQmean;
130 if (fRefPulserMasked) delete fRefPulserMasked;
131 if (fRefCETmean) delete fRefCETmean;
132 if (fRefCETrms) delete fRefCETrms;
133 if (fRefCEQmean) delete fRefCEQmean;
134 if (fRefCEMasked) delete fRefCEMasked;
135 if (fRefALTROFPED) delete fRefALTROFPED;
136 if (fRefALTROZsThr) delete fRefALTROZsThr;
137 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
138 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
139 if (fRefALTROMasked) delete fRefALTROMasked;
140 if (fRefCalibRaw) delete fRefCalibRaw;
141 if (fCurrentRefMap) delete fCurrentRefMap;
143 //_____________________________________________________________________________________
144 void AliTPCcalibDButil::UpdateFromCalibDB()
147 // Update pointers from calibDB
149 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
150 fCalibDB->UpdateNonRec(); // load all infromation now
151 fPadNoise=fCalibDB->GetPadNoise();
152 fPedestals=fCalibDB->GetPedestals();
153 fPulserTmean=fCalibDB->GetPulserTmean();
154 fPulserTrms=fCalibDB->GetPulserTrms();
155 fPulserQmean=fCalibDB->GetPulserQmean();
156 fCETmean=fCalibDB->GetCETmean();
157 fCETrms=fCalibDB->GetCETrms();
158 fCEQmean=fCalibDB->GetCEQmean();
159 fALTROMasked=fCalibDB->GetALTROMasked();
160 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
161 fCalibRaw=fCalibDB->GetCalibRaw();
162 fDataQA=fCalibDB->GetDataQA();
163 UpdatePulserOutlierMap();
164 // SetReferenceRun();
165 // UpdateRefDataFromOCDB();
167 //_____________________________________________________________________________________
168 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
169 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
172 // Process the CE data for this run
173 // the return TVectorD arrays contian the results of the fit
174 // noutliersCE contains the number of pads marked as outliers,
175 // not including masked and edge pads
178 //retrieve CE and ALTRO data
180 TString fitString(fitFormula);
181 fitString.ReplaceAll("++","#");
182 Int_t ndim=fitString.CountChar('#')+2;
183 fitResultsA.ResizeTo(ndim);
184 fitResultsC.ResizeTo(ndim);
193 if (outCE) out=outCE;
194 else out=new AliTPCCalPad("outCE","outCE");
195 AliTPCCalROC *rocMasked=0x0;
196 //loop over all channels
197 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
198 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
199 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
200 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
202 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
206 //add time offset to IROCs
207 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
208 rocData->Add(fIrocTimeOffset);
210 UInt_t nrows=rocData->GetNrows();
211 for (UInt_t irow=0;irow<nrows;++irow){
212 UInt_t npads=rocData->GetNPads(irow);
213 for (UInt_t ipad=0;ipad<npads;++ipad){
214 rocOut->SetValue(irow,ipad,0);
215 //exclude masked pads
216 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
217 rocOut->SetValue(irow,ipad,1);
220 //exclude first two rows in IROC and last two rows in OROC
222 if (irow<2) rocOut->SetValue(irow,ipad,1);
224 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
227 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
228 Float_t valTmean=rocData->GetValue(irow,ipad);
229 //exclude values that are exactly 0
230 if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
231 rocOut->SetValue(irow,ipad,1);
234 // exclude channels with too large variations
235 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
236 rocOut->SetValue(irow,ipad,1);
244 Float_t chi2Af,chi2Cf;
245 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
248 if (!outCE) delete out;
250 //_____________________________________________________________________________________
251 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
252 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
253 Float_t &driftTimeA, Float_t &driftTimeC )
256 // Calculate statistical information from the CE graphs for drift time and charge
260 vecTEntries.ResizeTo(72);
261 vecTMean.ResizeTo(72);
262 vecTRMS.ResizeTo(72);
263 vecTMedian.ResizeTo(72);
264 vecQEntries.ResizeTo(72);
265 vecQMean.ResizeTo(72);
266 vecQRMS.ResizeTo(72);
267 vecQMedian.ResizeTo(72);
278 TObjArray *arrT=fCalibDB->GetCErocTtime();
279 TObjArray *arrQ=fCalibDB->GetCErocQtime();
281 for (Int_t isec=0;isec<74;++isec){
282 TGraph *gr=(TGraph*)arrT->At(isec);
285 Int_t npoints = gr->GetN();
286 values.ResizeTo(npoints);
288 //skip first points, theres always some problems with finding the CE position
289 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
290 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
291 values[nused]=gr->GetY()[ipoint];
296 if (isec<72) vecTEntries[isec]= nused;
299 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
300 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
301 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
302 } else if (isec==72){
303 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
304 } else if (isec==73){
305 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
311 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
312 TGraph *gr=(TGraph*)arrQ->At(isec);
315 Int_t npoints = gr->GetN();
316 values.ResizeTo(npoints);
318 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
319 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
320 values[nused]=gr->GetY()[ipoint];
325 vecQEntries[isec]= nused;
327 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
328 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
329 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
335 //_____________________________________________________________________________________
336 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
337 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
338 Int_t &nonMaskedZero, Int_t &nNaN)
341 // process noise data
342 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
343 // OROCs small pads [2] and OROCs large pads [3]
344 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
345 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
348 //set proper size and reset
349 const UInt_t infoSize=4;
350 vNoiseMean.ResizeTo(infoSize);
351 vNoiseMeanSenRegions.ResizeTo(infoSize);
352 vNoiseRMS.ResizeTo(infoSize);
353 vNoiseRMSSenRegions.ResizeTo(infoSize);
355 vNoiseMeanSenRegions.Zero();
357 vNoiseRMSSenRegions.Zero();
361 TVectorD c(infoSize);
362 TVectorD cs(infoSize);
366 //retrieve noise and ALTRO data
367 if (!fPadNoise) return;
368 AliTPCCalROC *rocMasked=0x0;
369 //create IROC, OROC1, OROC2 and sensitive region masks
370 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
371 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
372 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
373 UInt_t nrows=noiseROC->GetNrows();
374 for (UInt_t irow=0;irow<nrows;++irow){
375 UInt_t npads=noiseROC->GetNPads(irow);
376 for (UInt_t ipad=0;ipad<npads;++ipad){
377 //don't use masked channels;
378 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
379 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
381 if (noiseVal<kAlmost0) {
386 if ( !(noiseVal<10000000) ){
387 // printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
391 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
392 Int_t masksen=1; // sensitive pards are not masked (0)
393 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
394 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
396 if (irow>19&&irow<46){
397 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
400 vNoiseMean[type]+=noiseVal;
401 vNoiseRMS[type]+=noiseVal*noiseVal;
404 vNoiseMeanSenRegions[type]+=noiseVal;
405 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
410 //define sensive regions
411 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
413 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
414 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
416 if ((Int_t)irow<par.GetNRowUp1()){
419 vNoiseMean[type]+=noiseVal;
420 vNoiseRMS[type]+=noiseVal*noiseVal;
423 vNoiseMeanSenRegions[type]+=noiseVal;
424 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
430 vNoiseMean[type]+=noiseVal;
431 vNoiseRMS[type]+=noiseVal*noiseVal;
434 vNoiseMeanSenRegions[type]+=noiseVal;
435 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
442 vNoiseMean[type]+=noiseVal;
443 vNoiseRMS[type]+=noiseVal*noiseVal;
446 vNoiseMeanSenRegions[type]+=noiseVal;
447 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
452 }//end loop sectors (rocs)
454 //calculate mean and RMS
455 const Double_t verySmall=0.0000000001;
456 for (UInt_t i=0;i<infoSize;++i){
463 // printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
464 mean=vNoiseMean[i]/c[i];
466 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
471 if (cs[i]>verySmall){
472 meanSen=vNoiseMeanSenRegions[i]/cs[i];
473 rmsSen=vNoiseRMSSenRegions[i];
474 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
476 vNoiseMeanSenRegions[i]=meanSen;
477 vNoiseRMSSenRegions[i]=rmsSen;
481 //_____________________________________________________________________________________
482 void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
488 // vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
492 const UInt_t infoSize = 72;
493 //reset counters to error number
494 vQaOcc.ResizeTo(infoSize);
496 vQaQtot.ResizeTo(infoSize);
498 vQaQmax.ResizeTo(infoSize);
501 //retrieve pulser and ALTRO data
505 AliInfo("No QA data");
508 if (fDataQA->GetEventCounter()<=0) {
510 AliInfo("No QA data");
511 return; // no data processed
516 TVectorD normOcc(infoSize);
517 TVectorD normQ(infoSize);
519 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
521 printf("Sector %d\n", isec);
523 AliTPCCalROC* occupancyROC = fDataQA->GetNoThreshold()->GetCalROC(isec);
524 AliTPCCalROC* nclusterROC = fDataQA->GetNLocalMaxima()->GetCalROC(isec);
525 AliTPCCalROC* qROC = fDataQA->GetMeanCharge()->GetCalROC(isec);
526 AliTPCCalROC* qmaxROC = fDataQA->GetMaxCharge()->GetCalROC(isec);
527 if (!occupancyROC) continue;
528 if (!nclusterROC) continue;
530 if (!qmaxROC) continue;
532 const UInt_t nchannels=occupancyROC->GetNchannels();
534 printf("Nchannels %d\n", nchannels);
536 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
538 vQaOcc[isec] += occupancyROC->GetValue(ichannel);
541 Float_t nClusters = nclusterROC->GetValue(ichannel);
542 normQ[isec] += nClusters;
543 vQaQtot[isec]+=nClusters*qROC->GetValue(ichannel);
544 vQaQmax[isec]+=nClusters*qmaxROC->GetValue(ichannel);
548 //calculate mean values
549 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
551 if (normOcc[isec]>0) vQaOcc[isec] /= normOcc[isec];
552 else vQaOcc[isec] = 0;
555 vQaQtot[isec] /= normQ[isec];
556 vQaQmax[isec] /= normQ[isec];
565 //_____________________________________________________________________________________
566 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
569 // Process the Pulser information
570 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
573 const UInt_t infoSize=4;
574 //reset counters to error number
575 vMeanTime.ResizeTo(infoSize);
578 TVectorD c(infoSize);
579 //retrieve pulser and ALTRO data
580 if (!fPulserTmean) return;
583 AliTPCCalROC *rocOut=0x0;
584 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
585 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
586 if (!tmeanROC) continue;
587 rocOut=fPulserOutlier->GetCalROC(isec);
588 UInt_t nchannels=tmeanROC->GetNchannels();
589 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
590 if (rocOut && rocOut->GetValue(ichannel)) continue;
591 Float_t val=tmeanROC->GetValue(ichannel);
593 vMeanTime[type]+=val;
598 for (UInt_t itype=0; itype<infoSize; ++itype){
599 if (c[itype]>0) vMeanTime[itype]/=c[itype];
600 else vMeanTime[itype]=0;
603 //_____________________________________________________________________________________
604 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
607 // Get Values from ALTRO configuration data
610 if (!fALTROMasked) return;
612 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
613 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
614 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
615 if (rocMasked->GetValue(ichannel)) ++nMasked;
619 //_____________________________________________________________________________________
620 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
623 // Proces Goofie values, return statistical information of the currently set goofieArray
624 // The meaning of the entries are given below
626 1 TPC_ANODE_I_A00_STAT
628 3 TPC_DVM_DriftVelocity
633 8 TPC_DVM_NumberOfSparks
634 9 TPC_DVM_PeakAreaFar
635 10 TPC_DVM_PeakAreaNear
636 11 TPC_DVM_PeakPosFar
637 12 TPC_DVM_PeakPosNear
643 18 TPC_DVM_TemperatureS1
647 vecEntries.ResizeTo(nsensors);
648 vecMedian.ResizeTo(nsensors);
649 vecMean.ResizeTo(nsensors);
650 vecRMS.ResizeTo(nsensors);
657 Double_t kEpsilon=0.0000000001;
658 Double_t kBig=100000000000.;
659 Int_t nsensors = fGoofieArray->NumSensors();
660 vecEntries.ResizeTo(nsensors);
661 vecMedian.ResizeTo(nsensors);
662 vecMean.ResizeTo(nsensors);
663 vecRMS.ResizeTo(nsensors);
665 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
666 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
667 if (gsensor && gsensor->GetGraph()){
668 Int_t npoints = gsensor->GetGraph()->GetN();
670 values.ResizeTo(npoints);
672 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
673 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
674 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
675 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
680 vecEntries[isensor]= nused;
682 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
683 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
684 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
689 //_____________________________________________________________________________________
690 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
693 // check the variations of the pedestal data to the reference pedestal data
694 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
697 TVectorF vThres(npar); //thresholds
698 Int_t nActive=0; //number of active channels
700 //reset and set thresholds
701 pedestalDeviations.ResizeTo(npar);
702 for (Int_t i=0;i<npar;++i){
703 pedestalDeviations.GetMatrixArray()[i]=0;
704 vThres.GetMatrixArray()[i]=(i+1)*.5;
706 //check all needed data is available
707 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
708 //loop over all channels
709 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
710 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
711 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
712 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
713 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
714 UInt_t nrows=mROC->GetNrows();
715 for (UInt_t irow=0;irow<nrows;++irow){
716 UInt_t npads=mROC->GetNPads(irow);
717 for (UInt_t ipad=0;ipad<npads;++ipad){
718 //don't use masked channels;
719 if (mROC ->GetValue(irow,ipad)) continue;
720 if (mRefROC->GetValue(irow,ipad)) continue;
721 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
722 for (Int_t i=0;i<npar;++i){
723 if (deviation>vThres[i])
724 ++pedestalDeviations.GetMatrixArray()[i];
731 for (Int_t i=0;i<npar;++i){
732 pedestalDeviations.GetMatrixArray()[i]/=nActive;
736 //_____________________________________________________________________________________
737 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
740 // check the variations of the noise data to the reference noise data
741 // thresholds are 5, 10, 15 and 20 percent respectively.
744 TVectorF vThres(npar); //thresholds
745 Int_t nActive=0; //number of active channels
747 //reset and set thresholds
748 noiseDeviations.ResizeTo(npar);
749 for (Int_t i=0;i<npar;++i){
750 noiseDeviations.GetMatrixArray()[i]=0;
751 vThres.GetMatrixArray()[i]=(i+1)*.05;
753 //check all needed data is available
754 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
755 //loop over all channels
756 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
757 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
758 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
759 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
760 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
761 UInt_t nrows=mROC->GetNrows();
762 for (UInt_t irow=0;irow<nrows;++irow){
763 UInt_t npads=mROC->GetNPads(irow);
764 for (UInt_t ipad=0;ipad<npads;++ipad){
765 //don't use masked channels;
766 if (mROC ->GetValue(irow,ipad)) continue;
767 if (mRefROC->GetValue(irow,ipad)) continue;
768 if (nRefROC->GetValue(irow,ipad)==0) continue;
769 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
770 for (Int_t i=0;i<npar;++i){
771 if (deviation>vThres[i])
772 ++noiseDeviations.GetMatrixArray()[i];
779 for (Int_t i=0;i<npar;++i){
780 noiseDeviations.GetMatrixArray()[i]/=nActive;
784 //_____________________________________________________________________________________
785 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
786 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
789 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
790 // thresholds are .5, 1, 5 and 10 percent respectively.
794 TVectorF vThres(npar); //thresholds
795 Int_t nActive=0; //number of active channels
797 //reset and set thresholds
798 pulserQdeviations.ResizeTo(npar);
799 for (Int_t i=0;i<npar;++i){
800 pulserQdeviations.GetMatrixArray()[i]=0;
805 vThres.GetMatrixArray()[0]=.005;
806 vThres.GetMatrixArray()[1]=.01;
807 vThres.GetMatrixArray()[2]=.05;
808 vThres.GetMatrixArray()[3]=.1;
809 //check all needed data is available
810 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
812 UpdateRefPulserOutlierMap();
813 //loop over all channels
814 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
815 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
816 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
817 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
818 // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
819 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
820 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
821 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
822 Float_t ptmean=ptROC->GetMean(oROC);
823 UInt_t nrows=mROC->GetNrows();
824 for (UInt_t irow=0;irow<nrows;++irow){
825 UInt_t npads=mROC->GetNPads(irow);
826 for (UInt_t ipad=0;ipad<npads;++ipad){
827 //don't use masked channels;
828 if (mROC ->GetValue(irow,ipad)) continue;
829 if (mRefROC->GetValue(irow,ipad)) continue;
830 //don't user edge pads
831 if (ipad==0||ipad==npads-1) continue;
833 Float_t pq=pqROC->GetValue(irow,ipad);
834 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
835 Float_t pt=ptROC->GetValue(irow,ipad);
836 // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
838 Float_t deviation=TMath::Abs(pq/pqRef-1);
839 for (Int_t i=0;i<npar;++i){
840 if (deviation>vThres[i])
841 ++pulserQdeviations.GetMatrixArray()[i];
843 if (pqRef>11&&pq<11) ++npadsOffAdd;
846 if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
852 for (Int_t i=0;i<npar;++i){
853 pulserQdeviations.GetMatrixArray()[i]/=nActive;
858 //_____________________________________________________________________________________
859 void AliTPCcalibDButil::UpdatePulserOutlierMap()
862 // Update the outlier map of the pulser data
864 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
866 //_____________________________________________________________________________________
867 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
870 // Update the outlier map of the pulser reference data
872 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
874 //_____________________________________________________________________________________
875 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
878 // Create a map that contains outliers from the Pulser calibration data.
879 // The outliers include masked channels, edge pads and pads with
880 // too large timing and charge variations.
881 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
882 // those do not contain masked and edge pads
886 pulOut->Multiply(0.);
890 AliTPCCalROC *rocMasked=0x0;
894 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
895 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
896 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
897 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
898 if (!tmeanROC||!qmeanROC) {
899 //reset outliers in this ROC
900 outROC->Multiply(0.);
903 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
905 // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
906 // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
907 UInt_t nrows=tmeanROC->GetNrows();
908 for (UInt_t irow=0;irow<nrows;++irow){
909 UInt_t npads=tmeanROC->GetNPads(irow);
910 for (UInt_t ipad=0;ipad<npads;++ipad){
911 Int_t outlier=0,masked=0;
912 Float_t q=qmeanROC->GetValue(irow,ipad);
913 Float_t t=tmeanROC->GetValue(irow,ipad);
914 //masked channels are outliers
915 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
916 //edge pads are outliers
917 if (ipad==0||ipad==npads-1) masked=1;
918 //channels with too large charge or timing deviation from the meadian are outliers
919 // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
920 if (q<fPulQminLimit && !masked) outlier=1;
922 if ( !(q<10000000) || !(t<10000000)) outlier=1;
923 outROC->SetValue(irow,ipad,outlier+masked);
924 fNpulserOutliers+=outlier;
929 //_____________________________________________________________________________________
930 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
933 // Create pad time0 object from pulser and/or CE data, depending on the selected model
934 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
935 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
936 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
938 // In case model 2 is invoked - gy arival time gradient is also returned
942 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
943 // decide between different models
944 if (model==0||model==1){
946 if (model==1) ProcessPulser(vMean);
947 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
948 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
949 if (!rocPulTmean) continue;
950 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
951 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
952 Float_t mean=rocPulTmean->GetMean(rocOut);
953 //treat case where a whole partition is masked
954 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
959 UInt_t nrows=rocTime0->GetNrows();
960 for (UInt_t irow=0;irow<nrows;++irow){
961 UInt_t npads=rocTime0->GetNPads(irow);
962 for (UInt_t ipad=0;ipad<npads;++ipad){
963 Float_t time=rocPulTmean->GetValue(irow,ipad);
964 //in case of an outlier pad use the mean of the altro values.
965 //This should be the most precise guess in that case.
966 if (rocOut->GetValue(irow,ipad)) {
967 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
968 if ( TMath::Abs(time)<kAlmost0 ) time=mean;
970 Float_t val=time-mean;
971 rocTime0->SetValue(irow,ipad,val);
975 } else if (model==2){
976 Double_t pgya,pgyc,pchi2a,pchi2c;
977 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
978 fCETmean->Add(padPulser,-1.);
980 AliTPCCalPad outCE("outCE","outCE");
982 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
983 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
984 // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
985 if (!padFit) { delete padPulser; return 0;}
988 fCETmean->Add(padPulser,1.);
989 padTime0->Add(fCETmean);
990 padTime0->Add(padFit,-1);
995 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
996 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
997 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
998 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
999 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
1000 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
1001 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
1002 Float_t mean=rocPulTmean->GetMean(rocOutPul);
1003 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
1004 UInt_t nrows=rocTime0->GetNrows();
1005 for (UInt_t irow=0;irow<nrows;++irow){
1006 UInt_t npads=rocTime0->GetNPads(irow);
1007 for (UInt_t ipad=0;ipad<npads;++ipad){
1008 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
1009 if (rocOutCE->GetValue(irow,ipad)){
1010 Float_t valOut=rocCEfit->GetValue(irow,ipad);
1011 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
1012 rocTime0->SetValue(irow,ipad,valOut);
1020 Double_t median = padTime0->GetMedian();
1021 padTime0->Add(-median); // normalize to median
1024 //_____________________________________________________________________________________
1025 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
1028 // GetMeanAlto information
1030 if (roc==0) return 0.;
1031 const Int_t sector=roc->GetSector();
1032 AliTPCROC *tpcRoc=AliTPCROC::Instance();
1033 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
1037 //loop over a small range around the requested pad (+-10 rows/pads)
1038 for (Int_t irow=row-10;irow<row+10;++irow){
1039 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
1040 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
1041 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
1042 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
1043 if (altroRoc!=altroCurr) continue;
1044 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
1045 Float_t val=roc->GetValue(irow,ipad);
1053 //_____________________________________________________________________________________
1054 void AliTPCcalibDButil::SetRefFile(const char* filename)
1057 // load cal pad objects form the reference file
1059 TDirectory *currDir=gDirectory;
1061 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
1062 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
1064 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
1065 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
1066 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
1068 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
1069 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
1070 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
1072 // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
1073 // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
1074 // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
1075 // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
1076 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
1080 //_____________________________________________________________________________________
1081 void AliTPCcalibDButil::UpdateRefDataFromOCDB()
1084 // set reference data from OCDB Reference map
1087 AliWarning("Referenc map not set!");
1092 AliCDBEntry* entry = 0x0;
1093 Bool_t hasAnyChanged=kFALSE;
1096 cdbPath="TPC/Calib/Pedestals";
1097 if (HasRefChanged(cdbPath.Data())){
1098 hasAnyChanged=kTRUE;
1099 //delete old entries
1100 if (fRefPedestals) delete fRefPedestals;
1101 if (fRefPedestalMasked) delete fRefPedestalMasked;
1102 fRefPedestals=fRefPedestalMasked=0x0;
1104 entry=GetRefEntry(cdbPath.Data());
1106 entry->SetOwner(kTRUE);
1107 fRefPedestals=GetRefCalPad(entry);
1109 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1114 cdbPath="TPC/Calib/PadNoise";
1115 if (HasRefChanged(cdbPath.Data())){
1116 hasAnyChanged=kTRUE;
1118 if (fRefPadNoise) delete fRefPadNoise;
1121 entry=GetRefEntry(cdbPath.Data());
1123 entry->SetOwner(kTRUE);
1124 fRefPadNoise=GetRefCalPad(entry);
1130 cdbPath="TPC/Calib/Pulser";
1131 if (HasRefChanged(cdbPath.Data())){
1132 hasAnyChanged=kTRUE;
1133 //delete old entries
1134 if (fRefPulserTmean) delete fRefPulserTmean;
1135 if (fRefPulserTrms) delete fRefPulserTrms;
1136 if (fRefPulserQmean) delete fRefPulserQmean;
1137 if (fRefPulserMasked) delete fRefPulserMasked;
1138 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1140 entry=GetRefEntry(cdbPath.Data());
1142 entry->SetOwner(kTRUE);
1143 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1144 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1145 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1147 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1152 cdbPath="TPC/Calib/CE";
1153 if (HasRefChanged(cdbPath.Data())){
1154 hasAnyChanged=kTRUE;
1155 //delete old entries
1156 if (fRefCETmean) delete fRefCETmean;
1157 if (fRefCETrms) delete fRefCETrms;
1158 if (fRefCEQmean) delete fRefCEQmean;
1159 if (fRefCEMasked) delete fRefCEMasked;
1160 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1162 entry=GetRefEntry(cdbPath.Data());
1164 entry->SetOwner(kTRUE);
1165 fRefCETmean=GetRefCalPad(entry,"CETmean");
1166 fRefCETrms=GetRefCalPad(entry,"CETrms");
1167 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1169 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1174 cdbPath="TPC/Calib/AltroConfig";
1175 if (HasRefChanged(cdbPath.Data())){
1176 hasAnyChanged=kTRUE;
1177 //delete old entries
1178 if (fRefALTROFPED) delete fRefALTROFPED;
1179 if (fRefALTROZsThr) delete fRefALTROZsThr;
1180 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1181 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1182 if (fRefALTROMasked) delete fRefALTROMasked;
1183 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1185 entry=GetRefEntry(cdbPath.Data());
1187 entry->SetOwner(kTRUE);
1188 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1189 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1190 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1191 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1192 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1199 cdbPath="TPC/Calib/Raw";
1200 if (HasRefChanged(cdbPath.Data())){
1201 hasAnyChanged=kTRUE;
1203 if (fRefCalibRaw) delete fRefCalibRaw;
1205 entry=GetRefEntry(cdbPath.Data());
1207 entry->SetOwner(kTRUE);
1208 TObjArray *arr=(TObjArray*)entry->GetObject();
1210 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1212 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1219 cdbPath="TPC/Calib/QA";
1220 if (HasRefChanged(cdbPath.Data())){
1221 hasAnyChanged=kTRUE;
1223 if (fRefDataQA) delete fRefDataQA;
1225 entry=GetRefEntry(cdbPath.Data());
1227 entry->SetOwner(kTRUE);
1228 fRefDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1230 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1232 fRefDataQA=(AliTPCdataQA*)fRefDataQA->Clone();
1239 //update current reference maps
1241 if (fCurrentRefMap) delete fCurrentRefMap;
1242 fCurrentRefMap=(TMap*)fRefMap->Clone();
1245 //_____________________________________________________________________________________
1246 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1249 // TObjArray object type case
1250 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1252 AliTPCCalPad *pad=0x0;
1253 TObjArray *arr=(TObjArray*)entry->GetObject();
1255 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1258 pad=(AliTPCCalPad*)arr->FindObject(objName);
1260 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
1263 return (AliTPCCalPad*)pad->Clone();
1265 //_____________________________________________________________________________________
1266 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1269 // AliTPCCalPad object type case
1270 // cast object to a calPad and store it in 'pad'
1272 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1274 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1277 pad=(AliTPCCalPad*)pad->Clone();
1280 //_____________________________________________________________________________________
1281 AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1284 // set altro masked channel map for 'cdbPath'
1286 AliTPCCalPad* pad=0x0;
1287 const Int_t run=GetReferenceRun(cdbPath);
1289 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1292 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1294 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1297 pad=GetRefCalPad(entry,"Masked");
1298 if (pad) pad->SetNameTitle(name,name);
1299 entry->SetOwner(kTRUE);
1303 //_____________________________________________________________________________________
1304 void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1306 // Get Reference map
1308 if (run<0) run=fCalibDB->GetRun();
1309 TString cdbPath="TPC/Calib/Ref";
1310 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1312 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
1316 entry->SetOwner(kTRUE);
1317 fRefMap=(TMap*)(entry->GetObject());
1318 AliCDBId &id=entry->GetId();
1319 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1321 //_____________________________________________________________________________________
1322 Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1325 // check whether a reference cdb entry has changed
1327 if (!fCurrentRefMap) return kTRUE;
1328 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1331 //_____________________________________________________________________________________
1332 AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1335 // get the reference AliCDBEntry for 'cdbPath'
1337 const Int_t run=GetReferenceRun(cdbPath);
1339 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1342 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1344 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1349 //_____________________________________________________________________________________
1350 Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
1352 // Get reference run number for the specified OCDB path
1354 if (!fCurrentRefMap) return -2;
1355 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1356 if (!str) return -2;
1357 return (Int_t)str->GetString().Atoi();
1359 //_____________________________________________________________________________________
1360 Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1362 // Get reference run number for the specified OCDB path
1364 if (!fRefMap) return -1;
1365 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1366 if (!str) return -1;
1367 return (Int_t)str->GetString().Atoi();
1369 //_____________________________________________________________________________________
1370 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
1372 // Author: marian.ivanov@cern.ch
1374 // Create outlier map for CE study
1376 // Return value - outlyer map
1377 // noutlyersCE - number of outlyers
1378 // minSignal - minimal total Q signal
1379 // cutRMSMin - minimal width of the signal in respect to the median
1380 // cutRMSMax - maximal width of the signal in respect to the median
1381 // cutMaxDistT - maximal deviation from time median per chamber
1383 // Outlyers criteria:
1384 // 0. Exclude masked pads
1385 // 1. Exclude first two rows in IROC and last two rows in OROC
1386 // 2. Exclude edge pads
1387 // 3. Exclude channels with too large variations
1388 // 4. Exclude pads with too small signal
1389 // 5. Exclude signal with outlyers RMS
1390 // 6. Exclude channels to far from the chamber median
1392 //create outlier map
1393 AliTPCCalPad *out=ceOut;
1394 if (!out) out= new AliTPCCalPad("outCE","outCE");
1395 AliTPCCalROC *rocMasked=0x0;
1396 if (!fCETmean) return 0;
1397 if (!fCETrms) return 0;
1398 if (!fCEQmean) return 0;
1400 //loop over all channels
1402 Double_t rmsMedian = fCETrms->GetMedian();
1403 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1404 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1405 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1406 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1407 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1408 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1409 Double_t trocMedian = rocData->GetMedian();
1411 if (!rocData || !rocCEQ || !rocCETrms) {
1412 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1418 UInt_t nrows=rocData->GetNrows();
1419 for (UInt_t irow=0;irow<nrows;++irow){
1420 UInt_t npads=rocData->GetNPads(irow);
1421 for (UInt_t ipad=0;ipad<npads;++ipad){
1422 rocOut->SetValue(irow,ipad,0);
1423 Float_t valTmean=rocData->GetValue(irow,ipad);
1424 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1425 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1426 //0. exclude masked pads
1427 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1428 rocOut->SetValue(irow,ipad,1);
1431 //1. exclude first two rows in IROC and last two rows in OROC
1433 if (irow<2) rocOut->SetValue(irow,ipad,1);
1435 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1437 //2. exclude edge pads
1438 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1439 //exclude values that are exactly 0
1440 if ( TMath::Abs(valTmean)<kAlmost0) {
1441 rocOut->SetValue(irow,ipad,1);
1444 //3. exclude channels with too large variations
1445 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1446 rocOut->SetValue(irow,ipad,1);
1450 //4. exclude channels with too small signal
1451 if (valQmean<minSignal) {
1452 rocOut->SetValue(irow,ipad,1);
1456 //5. exclude channels with too small rms
1457 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1458 rocOut->SetValue(irow,ipad,1);
1462 //6. exclude channels to far from the chamber median
1463 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1464 rocOut->SetValue(irow,ipad,1);
1475 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1477 // Author: marian.ivanov@cern.ch
1479 // Create outlier map for Pulser
1481 // Return value - outlyer map
1482 // noutlyersPulser - number of outlyers
1483 // cutTime - absolute cut - distance to the median of chamber
1484 // cutnRMSQ - nsigma cut from median q distribution per chamber
1485 // cutnRMSrms - nsigma cut from median rms distribution
1486 // Outlyers criteria:
1487 // 0. Exclude masked pads
1488 // 1. Exclude time outlyers (default 3 time bins)
1489 // 2. Exclude q outlyers (default 5 sigma)
1490 // 3. Exclude rms outlyers (default 5 sigma)
1492 AliTPCCalPad *out=pulserOut;
1493 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1494 AliTPCCalROC *rocMasked=0x0;
1495 if (!fPulserTmean) return 0;
1496 if (!fPulserTrms) return 0;
1497 if (!fPulserQmean) return 0;
1499 //loop over all channels
1501 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1502 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1503 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1504 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1505 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1506 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1508 Double_t rocMedianT = rocData->GetMedian();
1509 Double_t rocMedianQ = rocPulserQ->GetMedian();
1510 Double_t rocRMSQ = rocPulserQ->GetRMS();
1511 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1512 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1513 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1514 rocOut->SetValue(ichannel,0);
1515 Float_t valTmean=rocData->GetValue(ichannel);
1516 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1517 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1518 Float_t valMasked =0;
1519 if (rocMasked) valMasked = rocMasked->GetValue(ichannel);
1521 if (valMasked>0.5) isOut=1;
1522 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1523 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1524 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1525 rocOut->SetValue(ichannel,isOut);
1526 if (isOut) noutliersPulser++;
1533 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1535 // Author : Marian Ivanov
1536 // Create pad time0 correction map using information from the CE and from pulser
1539 // Return PadTime0 to be used for time0 relative alignment
1540 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1541 // using $ALICE_ROOT/TPC/script/gui application
1543 // fitResultsA - fitParameters A side
1544 // fitResultsC - fitParameters C side
1545 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1546 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1550 // 1. Find outlier map for CE
1551 // 2. Find outlier map for Pulser
1552 // 3. Replace outlier by median at given sector (median without outliers)
1553 // 4. Substract from the CE data pulser
1554 // 5. Fit the CE with formula
1555 // 5.1) (IROC-OROC) offset
1559 // 5.5) (IROC-OROC)*(lx-xmid)
1561 // 6. Substract gy fit dependence from the CE data
1562 // 7. Add pulser back to CE data
1563 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1564 // 9. return CE data
1566 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1567 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1570 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)";
1571 // output for fit formula
1572 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)";
1573 // gy part of formula
1574 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)";
1577 if (!fCETmean) return 0;
1578 Double_t pgya,pgyc,pchi2a,pchi2c;
1579 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1580 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1582 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1583 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1584 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1585 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1586 padPulser->SetName("padPulser");
1587 padPulserOut->SetName("padPulserOut");
1588 padCE->SetName("padCE");
1589 padCEIn->SetName("padCEIn");
1590 padCEOut->SetName("padCEOut");
1591 padOut->SetName("padOut");
1594 // make combined outlyers map
1595 // and replace outlyers in maps with median for chamber
1597 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1598 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1599 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1600 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1601 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1602 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1603 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1604 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1605 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1606 if (rocPulserOut->GetValue(ichannel)>0) {
1607 rocPulser->SetValue(ichannel,pulserMedian);
1608 rocOut->SetValue(ichannel,1);
1610 if (rocCEOut->GetValue(ichannel)>0) {
1611 rocCE->SetValue(ichannel,ceMedian);
1612 rocOut->SetValue(ichannel,1);
1617 // remove pulser time 0
1619 padCE->Add(padPulser,-1);
1624 Float_t chi2Af,chi2Cf;
1625 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1629 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1630 padCEFitGY->SetName("padCEFitGy");
1632 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1633 padCEFit->SetName("padCEFit");
1635 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1636 padCEDiff->SetName("padCEDiff");
1637 padCEDiff->Add(padCEFit,-1.);
1640 padCE->Add(padCEFitGY,-1.);
1642 padCE->Add(padPulser,1.);
1643 Double_t padmedian = padCE->GetMedian();
1644 padCE->Add(-padmedian); // normalize to median
1646 // Replace outliers by fit value - median of diff per given chamber -GY fit
1648 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1649 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1650 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1651 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1652 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1653 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1655 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1656 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1657 if (rocOut->GetValue(ichannel)==0) continue;
1658 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1659 rocCE->SetValue(ichannel,value);
1665 //dump to the file - result can be visualized
1666 AliTPCPreprocessorOnline preprocesor;
1667 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1668 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1669 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1670 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1672 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1673 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1675 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1676 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1677 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1678 preprocesor.DumpToFile(dumpfile);
1681 delete padPulserOut;
1694 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1696 // find the closest point to xref in x direction
1697 // return dx and value
1701 if(!graph) return 0;
1702 if(graph->GetN() < 1) return 0;
1705 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1706 if (index<0) index=0;
1707 if(graph->GetN()==1) {
1708 dx = xref-graph->GetX()[index];
1711 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1712 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1713 dx = xref-graph->GetX()[index];
1715 y = graph->GetY()[index];
1719 Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1721 // Get the correction of the trigger offset
1722 // combining information from the laser track calibration
1723 // and from cosmic calibration
1726 // timeStamp - tim stamp in seconds
1727 // deltaT - integration period to calculate offset
1728 // deltaTLaser -max validity of laser data
1729 // valType - 0 - median, 1- mean
1731 // Integration vaues are just recomendation - if not possible to get points
1732 // automatically increase the validity by factor 2
1733 // (recursive algorithm until one month of data taking)
1736 const Float_t kLaserCut=0.0005;
1737 const Int_t kMaxPeriod=3600*24*30*12; // one year max
1738 const Int_t kMinPoints=20;
1740 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1742 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1744 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1745 if (!array) return 0;
1747 TGraphErrors *laserA[3]={0,0,0};
1748 TGraphErrors *laserC[3]={0,0,0};
1749 TGraphErrors *cosmicAll=0;
1750 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1751 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1752 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1755 if (!cosmicAll) return 0;
1756 Int_t nmeasC=cosmicAll->GetN();
1757 Float_t *tdelta = new Float_t[nmeasC];
1759 for (Int_t i=0;i<nmeasC;i++){
1760 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1761 Float_t ccosmic=cosmicAll->GetY()[i];
1762 Double_t yA=0,yC=0,dA=0,dC=0;
1763 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1764 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1765 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1766 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1768 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1770 if (TMath::Abs(yA-yC)<kLaserCut) {
1773 if (i%2==0) claser=yA;
1774 if (i%2==1) claser=yC;
1776 tdelta[nused]=ccosmic-claser;
1779 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1780 if (nused<kMinPoints) {
1782 printf("AliFatal: No time offset calibration available\n");
1785 Double_t median = TMath::Median(nused,tdelta);
1786 Double_t mean = TMath::Mean(nused,tdelta);
1788 return (valType==0) ? median:mean;
1791 Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1793 // Get the correction of the drift velocity
1794 // combining information from the laser track calibration
1795 // and from cosmic calibration
1797 // dist - return value - distance to closest point in graph
1799 // timeStamp - tim stamp in seconds
1800 // deltaT - integration period to calculate time0 offset
1801 // deltaTLaser -max validity of laser data
1802 // valType - 0 - median, 1- mean
1804 // Integration vaues are just recomendation - if not possible to get points
1805 // automatically increase the validity by factor 2
1806 // (recursive algorithm until one month of data taking)
1810 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1812 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1814 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1815 if (!array) return 0;
1816 TGraphErrors *cosmicAll=0;
1817 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1818 if (!cosmicAll) return 0;
1820 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1822 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1823 Double_t vcosmic = AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1824 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1825 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1832 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1833 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1834 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1836 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1837 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1838 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1839 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1841 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1842 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1848 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1851 // Create a default name for the gui file
1854 return Form("guiRefTreeRun%s.root",GetRefValidity());
1857 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1860 // Create a gui reference tree
1861 // if dirname and filename are empty default values will be used
1862 // this is the recommended way of using this function
1863 // it allows to check whether a file with the given run validity alredy exists
1865 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1866 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1870 TString file=filename;
1871 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1873 AliTPCPreprocessorOnline prep;
1874 //noise and pedestals
1875 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1876 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1877 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1879 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1880 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1881 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1882 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1884 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1885 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1886 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1887 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1889 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1890 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1891 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1892 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1893 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1895 AliTPCdataQA *dataQA=fRefDataQA;
1897 if (dataQA->GetNLocalMaxima())
1898 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1899 if (dataQA->GetMaxCharge())
1900 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1901 if (dataQA->GetMeanCharge())
1902 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1903 if (dataQA->GetNoThreshold())
1904 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1905 if (dataQA->GetNTimeBins())
1906 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1907 if (dataQA->GetNPads())
1908 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1909 if (dataQA->GetTimePosition())
1910 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1912 prep.DumpToFile(file.Data());
1916 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1918 // Get the correction of the drift velocity using the offline laser tracks calbration
1921 // timeStamp - tim stamp in seconds
1922 // deltaT - integration period to calculate time0 offset
1923 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1924 // Note in case no data form both A and C side - the value from active side used
1925 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1927 return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1930 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1932 // Get the correction of the drift velocity using the online laser tracks calbration
1935 // timeStamp - tim stamp in seconds
1936 // deltaT - integration period to calculate time0 offset
1937 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1938 // Note in case no data form both A and C side - the value from active side used
1939 TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1941 Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1942 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1943 if (!param) return 0;
1945 //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1946 dv*=param->GetDriftV()/2.61301900000000000e+06;
1947 if (dv>1e-20) dv=1/dv-1;
1950 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1952 AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
1953 AliDCSSensor *press = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
1959 AliTPCCalibVdrift corr(temp,press,0);
1960 corrPTA=corr.GetPTRelative(timeStamp,0);
1961 corrPTC=corr.GetPTRelative(timeStamp,1);
1964 if (side==0) dv -= corrPTA;
1965 if (side==1) dv -= corrPTC;
1966 if (side==2) dv -= (corrPTA+corrPTC)/2;
1971 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1972 Int_t side, TObjArray * const array){
1974 // common drift velocity retrieval for online and offline method
1976 TGraphErrors *grlaserA=0;
1977 TGraphErrors *grlaserC=0;
1978 Double_t vlaserA=0, vlaserC=0;
1979 if (!array) return 0;
1980 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1981 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1983 if (grlaserA && grlaserA->GetN()>0) {
1984 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1985 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1986 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1988 if (grlaserC && grlaserC->GetN()>0) {
1989 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1990 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1991 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1993 if (side==0) return vlaserA;
1994 if (side==1) return vlaserC;
1995 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1996 if (!grlaserA) return vlaserC;
1997 if (!grlaserC) return vlaserA;
2002 Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
2004 // Get the correction of the drift velocity using the CE laser data
2005 // combining information from the CE, laser track calibration
2006 // and P/T calibration
2009 // timeStamp - tim stamp in seconds
2010 // deltaT - integration period to calculate time0 offset
2011 // side - 0 - A side, 1 - C side, 2 - mean from both sides
2012 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
2013 if (!arrT) return 0;
2014 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
2015 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2016 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2019 Double_t corrPTA = 0, corrPTC=0;
2020 Double_t ltime0A = 0, ltime0C=0;
2022 Double_t corrA=0, corrC=0;
2023 Double_t timeA=0, timeC=0;
2024 const Double_t kEpsilon = 0.00001;
2025 TGraph *graphA = (TGraph*)arrT->At(72);
2026 TGraph *graphC = (TGraph*)arrT->At(73);
2027 if (!graphA && !graphC) return 0.;
2028 if (graphA &&graphA->GetN()>0) {
2029 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2030 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2031 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2032 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2033 if(ltime0A < kEpsilon) return 0;
2034 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
2035 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2038 if (graphC&&graphC->GetN()>0){
2039 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2040 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2041 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2042 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2043 if(ltime0C < kEpsilon) return 0;
2044 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
2045 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2049 if (side ==0 ) return corrA;
2050 if (side ==1 ) return corrC;
2051 Double_t corrM= (corrA+corrC)*0.5;
2052 if (!graphA) corrM=corrC;
2053 if (!graphC) corrM=corrA;
2057 Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2059 // return drift velocity using the TPC-ITS matchin method
2060 // return also distance to the closest point
2062 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2063 TGraphErrors *graph=0;
2065 if (!array) return 0;
2067 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2068 if (!graph) return 0;
2070 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2071 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2075 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2077 // Get time dependent time 0 (trigger delay in cm) correction
2079 // timestamp - timestamp
2082 // Notice - Extrapolation outside of calibration range - using constant function
2084 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2085 TGraphErrors *graph=0;
2087 if (!array) return 0;
2088 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
2089 if (!graph) return 0;
2091 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2092 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2100 Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2102 // VERY obscure method - we need something in framework
2103 // Find the TPC runs with temperature OCDB entry
2104 // cache the start and end of the run
2106 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2107 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2108 if (!storage) return 0;
2109 TString path=storage->GetURI();
2113 if (path.Contains("local")){ // find the list if local system
2114 path.ReplaceAll("local://","");
2115 path+="TPC/Calib/Temperature";
2116 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
2118 runsT=gSystem->GetFromPipe(command);
2120 TObjArray *arr= runsT.Tokenize("r");
2123 TArrayI indexes(arr->GetEntries());
2124 TArrayI runs(arr->GetEntries());
2126 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2127 Int_t irunN = atoi(arr->At(irun)->GetName());
2128 if (irunN<startRun) continue;
2129 if (irunN>stopRun) continue;
2130 runs[naccept]=irunN;
2134 fRunsStart.Set(fRuns.fN);
2135 fRunsStop.Set(fRuns.fN);
2136 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2137 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
2140 AliCDBEntry * entry = 0;
2141 {for (Int_t irun=0;irun<fRuns.fN; irun++){
2142 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2143 if (!entry) continue;
2144 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2145 if (!tmpRun) continue;
2146 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2147 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2148 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2154 Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2156 // binary search - find the run for given time stamp
2158 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2159 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2161 for (Int_t index=index0; index<=index1; index++){
2162 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2164 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2167 if (cindex<0) cindex =(index0+index1)/2;
2171 return fRuns[cindex];
2178 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2180 // filter outlyer measurement
2181 // Only points around median +- sigmaCut filtered
2182 if (!graph) return 0;
2184 Int_t npoints0 = graph->GetN();
2189 if (npoints0<kMinPoints) return 0;
2191 Double_t *outx=new Double_t[npoints0];
2192 Double_t *outy=new Double_t[npoints0];
2193 for (Int_t iter=0; iter<3; iter++){
2195 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2196 if (graph->GetY()[ipoint]==0) continue;
2197 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2198 outx[npoints] = graph->GetX()[ipoint];
2199 outy[npoints] = graph->GetY()[ipoint];
2202 if (npoints<=1) break;
2203 medianY =TMath::Median(npoints,outy);
2204 rmsY =TMath::RMS(npoints,outy);
2207 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2214 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2216 // filter outlyer measurement
2217 // Only points around median +- cut filtered
2218 if (!graph) return 0;
2220 Int_t npoints0 = graph->GetN();
2225 if (npoints0<kMinPoints) return 0;
2227 Double_t *outx=new Double_t[npoints0];
2228 Double_t *outy=new Double_t[npoints0];
2229 for (Int_t iter=0; iter<3; iter++){
2231 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2232 if (graph->GetY()[ipoint]==0) continue;
2233 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2234 outx[npoints] = graph->GetX()[ipoint];
2235 outy[npoints] = graph->GetY()[ipoint];
2238 if (npoints<=1) break;
2239 medianY =TMath::Median(npoints,outy);
2240 rmsY =TMath::RMS(npoints,outy);
2243 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2251 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2253 // filter outlyer measurement
2254 // Only points with normalized errors median +- sigmaCut filtered
2256 Int_t kMinPoints=10;
2257 Int_t npoints0 = graph->GetN();
2259 Float_t medianErr=0, rmsErr=0;
2262 if (npoints0<kMinPoints) return 0;
2264 Double_t *outx=new Double_t[npoints0];
2265 Double_t *outy=new Double_t[npoints0];
2266 Double_t *erry=new Double_t[npoints0];
2267 Double_t *nerry=new Double_t[npoints0];
2268 Double_t *errx=new Double_t[npoints0];
2270 for (Int_t iter=0; iter<3; iter++){
2272 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2273 nerry[npoints] = graph->GetErrorY(ipoint);
2274 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2275 erry[npoints] = graph->GetErrorY(ipoint);
2276 outx[npoints] = graph->GetX()[ipoint];
2277 outy[npoints] = graph->GetY()[ipoint];
2278 errx[npoints] = graph->GetErrorY(ipoint);
2281 if (npoints==0) break;
2282 medianErr=TMath::Median(npoints,erry);
2283 medianY =TMath::Median(npoints,outy);
2284 rmsErr =TMath::RMS(npoints,erry);
2286 TGraphErrors *graphOut=0;
2287 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2297 void AliTPCcalibDButil::Sort(TGraph *graph){
2299 // sort array - neccessay for approx
2301 Int_t npoints = graph->GetN();
2302 Int_t *indexes=new Int_t[npoints];
2303 Double_t *outx=new Double_t[npoints];
2304 Double_t *outy=new Double_t[npoints];
2305 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2306 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2307 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2308 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2309 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2315 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2317 // smmoth graph - mean on the interval
2320 Int_t npoints = graph->GetN();
2321 Double_t *outy=new Double_t[npoints];
2323 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2324 Double_t lx=graph->GetX()[ipoint];
2325 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2326 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2327 if (index0<0) index0=0;
2328 if (index1>=npoints-1) index1=npoints-1;
2329 if ((index1-index0)>1){
2330 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2332 outy[ipoint]=graph->GetY()[ipoint];
2335 // TLinearFitter fitter(3,"pol2");
2336 // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2337 // Double_t lx=graph->GetX()[ipoint];
2338 // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2339 // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2340 // if (index0<0) index0=0;
2341 // if (index1>=npoints-1) index1=npoints-1;
2342 // fitter.ClearPoints();
2343 // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2344 // if ((index1-index0)>1){
2345 // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2347 // outy[ipoint]=graph->GetY()[ipoint];
2353 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2354 graph->GetY()[ipoint] = outy[ipoint];
2359 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2361 // Use constant interpolation outside of range
2364 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2368 if (graph->GetN()<1){
2369 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
2374 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2375 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2377 // printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
2379 if(graph->GetN()==1)
2380 return graph->Eval(graph->GetX()[0]);
2383 return graph->Eval(xref);
2386 Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2388 // Use constant interpolation outside of range also for spline fits
2391 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2394 if (graph->GetKnots()<1){
2395 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2398 if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2399 if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2400 return graph->Eval( xref);
2403 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2405 // Filter DCS sensor information
2406 // ymin - minimal value
2408 // maxdy - maximal deirivative
2409 // sigmaCut - cut on values and derivative in terms of RMS distribution
2410 // Return value - accepted fraction
2414 // 0. Calculate median and rms of values in specified range
2415 // 1. Filter out outliers - median+-sigmaCut*rms
2416 // values replaced by median
2418 AliSplineFit * fit = sensor->GetFit();
2419 if (!fit) return 0.;
2420 Int_t nknots = fit->GetKnots();
2427 Double_t *yin0 = new Double_t[nknots];
2428 Double_t *yin1 = new Double_t[nknots];
2431 for (Int_t iknot=0; iknot< nknots; iknot++){
2432 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2433 yin0[naccept] = fit->GetY0()[iknot];
2434 yin1[naccept] = fit->GetY1()[iknot];
2435 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2447 Double_t medianY0=0, medianY1=0;
2448 Double_t rmsY0 =0, rmsY1=0;
2449 medianY0 = TMath::Median(naccept, yin0);
2450 medianY1 = TMath::Median(naccept, yin1);
2451 rmsY0 = TMath::RMS(naccept, yin0);
2452 rmsY1 = TMath::RMS(naccept, yin1);
2455 // 1. Filter out outliers - median+-sigmaCut*rms
2456 // values replaced by median
2457 // if replaced the derivative set to 0
2459 for (Int_t iknot=0; iknot< nknots; iknot++){
2461 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2462 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2463 if (nknots<2) fit->GetY1()[iknot]=0;
2464 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2466 fit->GetY0()[iknot]=medianY0;
2467 fit->GetY1()[iknot]=0;
2474 return Float_t(naccept)/Float_t(nknots);
2477 Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2479 // Filter temperature array
2480 // tempArray - array of temperatures -
2481 // ymin - minimal accepted temperature - default 15
2482 // ymax - maximal accepted temperature - default 22
2483 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2484 // return value - fraction of filtered sensors
2485 const Double_t kMaxDy=0.1;
2486 Int_t nsensors=tempArray->NumSensors();
2487 if (nsensors==0) return 0.;
2489 for (Int_t isensor=0; isensor<nsensors; isensor++){
2490 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2491 if (!sensor) continue;
2492 //printf("%d\n",isensor);
2493 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2494 if (sensor->GetFit()==0){
2496 tempArray->RemoveSensorNum(isensor);
2501 return Float_t(naccept)/Float_t(nsensors);
2505 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2508 // Input parameters:
2509 // deltaT - smoothing window (in seconds)
2510 // cutAbs - max distance of the time info to the median (in time bins)
2511 // cutSigma - max distance (in the RMS)
2512 // pcstream - optional debug streamer to store original and filtered info
2513 // Hardwired parameters:
2514 // kMinPoints =10; // minimal number of points to define the CE
2515 // kMinSectors=12; // minimal number of sectors to define sideCE
2517 // 0. Filter almost emty graphs (kMinPoints=10)
2518 // 1. calculate median and RMS per side
2519 // 2. Filter graphs - in respect with side medians
2520 // - cutAbs and cutDelta used
2521 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2522 // 4. Calculate mean for A side and C side
2524 const Int_t kMinPoints =10; // minimal number of points to define the CE
2525 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2526 const Int_t kMinTime =400; // minimal arrival time of CE
2527 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2529 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2530 if (!cearray) return;
2535 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2536 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2537 if ( tempMapCE && cavernPressureCE){
2539 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2540 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2541 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2544 // recalculate P/T correction map for time of the CE
2545 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2546 driftCalib->SetName("driftPTCE");
2547 driftCalib->SetTitle("driftPTCE");
2548 cearray->AddLast(driftCalib);
2552 // 0. Filter almost emty graphs
2555 for (Int_t i=0; i<72;i++){
2556 TGraph *graph= (TGraph*)arrT->At(i);
2557 if (!graph) continue;
2559 if (graph->GetN()<kMinPoints){
2561 delete graph; // delete empty graph
2564 if (tmin<0) tmin = graph->GetX()[0];
2565 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2567 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2568 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2571 // 1. calculate median and RMS per side
2573 TArrayF arrA(100000), arrC(100000);
2575 Double_t medianA=0, medianC=0;
2576 Double_t rmsA=0, rmsC=0;
2577 for (Int_t isec=0; isec<72;isec++){
2578 TGraph *graph= (TGraph*)arrT->At(isec);
2579 if (!graph) continue;
2580 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2581 if (graph->GetY()[ipoint]<kMinTime) continue;
2582 if (nA>=arrA.fN) arrA.Set(nA*2);
2583 if (nC>=arrC.fN) arrC.Set(nC*2);
2584 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2585 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2589 medianA=TMath::Median(nA,arrA.fArray);
2590 rmsA =TMath::RMS(nA,arrA.fArray);
2593 medianC=TMath::Median(nC,arrC.fArray);
2594 rmsC =TMath::RMS(nC,arrC.fArray);
2597 // 2. Filter graphs - in respect with side medians
2599 TArrayD vecX(100000), vecY(100000);
2600 for (Int_t isec=0; isec<72;isec++){
2601 TGraph *graph= (TGraph*)arrT->At(isec);
2602 if (!graph) continue;
2603 Double_t median = (isec%36<18) ? medianA: medianC;
2604 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2606 // for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2607 for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
2608 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2609 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2610 vecX[naccept]= graph->GetX()[ipoint];
2611 vecY[naccept]= graph->GetY()[ipoint];
2614 if (naccept<kMinPoints){
2615 arrT->AddAt(0,isec);
2616 delete graph; // delete empty graph
2619 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2621 arrT->AddAt(graph2,isec);
2624 // 3. Cut in respect wit the graph median
2626 for (Int_t i=0; i<72;i++){
2627 TGraph *graph= (TGraph*)arrT->At(i);
2628 if (!graph) continue;
2632 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2633 if (!graphTS0) continue;
2634 if (graphTS0->GetN()<kMinPoints) {
2640 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2641 if (!graphTS) continue;
2643 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2645 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2646 (*pcstream)<<"filterCE"<<
2651 "graphTS0.="<<graphTS0<<
2652 "graphTS.="<<graphTS<<
2656 arrT->AddAt(graphTS,i);
2660 // Recalculate the mean time A side C side
2662 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2663 Int_t meanPoints=(nA+nC)/72; // mean number of points
2664 for (Int_t itime=0; itime<200; itime++){
2666 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2667 for (Int_t i=0; i<72;i++){
2668 TGraph *graph= (TGraph*)arrT->At(i);
2669 if (!graph) continue;
2670 if (graph->GetN()<(meanPoints/4)) continue;
2671 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2672 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2676 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2677 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2678 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2679 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2682 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2683 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2685 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2686 (*pcstream)<<"filterAC"<<
2695 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2696 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2697 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2698 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2699 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2700 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2703 if (nA<kMinSectors) arrT->AddAt(0,72);
2704 else arrT->AddAt(graphTSA,72);
2705 if (nC<kMinSectors) arrT->AddAt(0,73);
2706 else arrT->AddAt(graphTSC,73);
2710 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2712 // Filter Drift velocity measurement using the tracks
2713 // 0. remove outlyers - error based
2717 const Int_t kMinPoints=1; // minimal number of points to define value
2718 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2721 for (Int_t i=0; i<arrT->GetEntries();i++){
2722 TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
2723 if (!graph) continue;
2724 if (graph->GetN()<kMinPoints){
2729 TGraphErrors *graph2 = NULL;
2730 if(graph->GetN()<10) {
2731 graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2733 delete graph; arrT->AddAt(0,i); continue;
2737 graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2739 delete graph; arrT->AddAt(0,i); continue;
2742 if (graph2->GetN()<1) {
2743 delete graph; arrT->AddAt(0,i); continue;
2745 graph2->SetName(graph->GetName());
2746 graph2->SetTitle(graph->GetTitle());
2747 arrT->AddAt(graph2,i);
2749 (*pcstream)<<"filterTracks"<<
2754 "graph2.="<<graph2<<
2765 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2768 // get laser time offset
2769 // median around timeStamp+-deltaT
2770 // QA - chi2 needed for later usage - to be added
2771 // - currently cut on error
2774 Double_t kMinDelay=0.01;
2775 Double_t kMinDelayErr=0.0001;
2777 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2778 if (!array) return 0;
2779 TGraphErrors *tlaser=0;
2781 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2782 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2784 if (!tlaser) return 0;
2785 Int_t npoints0= tlaser->GetN();
2786 if (npoints0==0) return 0;
2787 Double_t *xlaser = new Double_t[npoints0];
2788 Double_t *ylaser = new Double_t[npoints0];
2790 for (Int_t i=0;i<npoints0;i++){
2792 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2793 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2794 xlaser[npoints]=tlaser->GetX()[npoints];
2795 ylaser[npoints]=tlaser->GetY()[npoints];
2800 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2801 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2802 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2803 if (index0<0) index0=0;
2804 if (index1>=npoints-1) index1=npoints-1;
2805 if (index1-index0<kMinPoints) {
2811 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2812 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2821 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2823 // Filter Goofie data
2824 // goofieArray - points will be filtered
2825 // deltaT - smmothing time window
2826 // cutSigma - outler sigma cut in rms
2827 // minVn, maxVd- range absolute cut for variable vd/pt
2830 // Ignore goofie if not enough points
2832 const Int_t kMinPoints = 3;
2835 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2836 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2837 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2838 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2839 if (!graphvd) return;
2840 if (graphvd->GetN()<kMinPoints){
2842 goofieArray->GetSensorNum(2)->SetGraph(0);
2846 // 1. Caluclate medians of critical variables
2852 Double_t medianpt=0;
2853 Double_t medianvd=0, sigmavd=0;
2854 Double_t medianan=0;
2855 Double_t medianaf=0;
2856 Int_t entries=graphvd->GetN();
2857 Double_t yvdn[10000];
2860 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2861 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2862 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2863 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2864 yvdn[nvd++]=graphvd->GetY()[ipoint];
2866 if (nvd<kMinPoints){
2868 goofieArray->GetSensorNum(2)->SetGraph(0);
2872 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2873 if (nuni>=kMinPoints){
2874 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2876 medianvd = TMath::Median(nvd, yvdn);
2879 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2880 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2881 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2882 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2883 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2884 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2892 // 2. Make outlyer graph
2895 TGraph graphOut(*graphvd);
2896 for (Int_t i=0; i<entries;i++){
2898 Bool_t isOut=kFALSE;
2899 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2900 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2902 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2904 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2905 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2906 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2907 graphOut.GetY()[i]= (isOut)?1:0;
2910 if (nOK<kMinPoints) {
2912 goofieArray->GetSensorNum(2)->SetGraph(0);
2916 // 3. Filter out outlyers - and smooth
2918 TVectorF vmedianArray(goofieArray->NumSensors());
2919 TVectorF vrmsArray(goofieArray->NumSensors());
2920 Double_t xnew[10000];
2921 Double_t ynew[10000];
2923 junk.SetOwner(kTRUE);
2927 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2929 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2930 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2932 if (!sensor) continue;
2933 graphOld = sensor->GetGraph();
2935 sensor->SetGraph(0);
2937 for (Int_t i=0;i<entries;i++){
2938 if (graphOut.GetY()[i]>0.5) continue;
2939 xnew[nused]=graphOld->GetX()[i];
2940 ynew[nused]=graphOld->GetY()[i];
2943 graphNew = new TGraph(nused,xnew,ynew);
2944 junk.AddLast(graphNew);
2945 junk.AddLast(graphOld);
2947 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2949 junk.AddLast(graphNew0);
2950 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2952 junk.AddLast(graphNew1);
2953 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2955 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2956 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2957 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2958 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2959 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2960 vmedianArray[isensor]=median;
2966 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2967 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2968 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2969 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2970 (*pcstream)<<"goofieA"<<
2971 Form("isOK_%d.=",isensor)<<isOK<<
2972 Form("s_%d.=",isensor)<<sensor<<
2973 Form("gr_%d.=",isensor)<<graphOld<<
2974 Form("gr0_%d.=",isensor)<<graphNew0<<
2975 Form("gr1_%d.=",isensor)<<graphNew1<<
2976 Form("gr2_%d.=",isensor)<<graphNew2;
2977 if (isOK) sensor->SetGraph(graphNew2);
2979 (*pcstream)<<"goofieA"<<
2980 "vmed.="<<&vmedianArray<<
2981 "vrms.="<<&vrmsArray<<
2983 junk.Delete(); // delete temoprary graphs
2991 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2993 // Make a statistic matrix
2994 // Input parameters:
2995 // array - TObjArray of AliRelKalmanAlign
2996 // minFraction - minimal ration of accepted tracks
2997 // minStat - minimal statistic (number of accepted tracks)
2998 // maxvd - maximal deviation for the 1
3000 // columns - Mean, Median, RMS
3001 // row - parameter type (rotation[3], translation[3], drift[3])
3002 if (!array) return 0;
3003 if (array->GetEntries()<=0) return 0;
3004 // Int_t entries = array->GetEntries();
3005 Int_t entriesFast = array->GetEntriesFast();
3007 TVectorD *valArray[9];
3008 for (Int_t i=0; i<9; i++){
3009 valArray[i] = new TVectorD(entriesFast);
3012 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3013 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3014 if (!kalman) continue;
3015 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3016 if (kalman->GetNUpdates()<minStat) continue;
3017 if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
3018 kalman->GetState(state);
3019 for (Int_t ipar=0; ipar<9; ipar++)
3020 (*valArray[ipar])[naccept]=state[ipar];
3023 //if (naccept<2) return 0;
3024 if (naccept<1) return 0;
3025 TMatrixD *pstat=new TMatrixD(9,3);
3026 TMatrixD &stat=*pstat;
3027 for (Int_t ipar=0; ipar<9; ipar++){
3028 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3029 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3030 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3036 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
3038 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3040 // array - input array
3041 // stat - mean parameters statistic
3043 // sigmaCut - maximal allowed deviation from mean in terms of RMS
3044 if (!array) return 0;
3045 if (array->GetEntries()<=0) return 0;
3046 if (!(&stat)) return 0;
3047 // error increase in 1 hour
3048 const Double_t kerrsTime[9]={
3049 0.00001, 0.00001, 0.00001,
3050 0.001, 0.001, 0.001,
3051 0.002, 0.01, 0.001};
3054 Int_t entries = array->GetEntriesFast();
3055 TObjArray *sArray= new TObjArray(entries);
3056 AliRelAlignerKalman * sKalman =0;
3058 for (Int_t i=0; i<entries; i++){
3059 Int_t index=(direction)? entries-i-1:i;
3060 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3061 if (!kalman) continue;
3063 kalman->GetState(state);
3064 for (Int_t ipar=0; ipar<9; ipar++){
3065 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3067 if (!sKalman &&isOK) {
3068 sKalman=new AliRelAlignerKalman(*kalman);
3069 sKalman->SetRejectOutliers(kFALSE);
3070 sKalman->SetRunNumber(kalman->GetRunNumber());
3071 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3073 if (!sKalman) continue;
3074 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
3075 for (Int_t ipar=0; ipar<9; ipar++){
3076 // (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3077 // (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3078 // (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3079 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3081 sKalman->SetRunNumber(kalman->GetRunNumber());
3082 if (!isOK) sKalman->SetRunNumber(0);
3083 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3084 if (!isOK) continue;
3085 sKalman->SetRejectOutliers(kFALSE);
3086 sKalman->SetRunNumber(kalman->GetRunNumber());
3087 sKalman->SetTimeStamp(kalman->GetTimeStamp());
3088 sKalman->Merge(kalman);
3089 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3095 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
3097 // Merge 2 RelKalman arrays
3099 // arrayP - rel kalman in direction plus
3100 // arrayM - rel kalman in direction minus
3101 if (!arrayP) return 0;
3102 if (arrayP->GetEntries()<=0) return 0;
3103 if (!arrayM) return 0;
3104 if (arrayM->GetEntries()<=0) return 0;
3106 Int_t entries = arrayP->GetEntriesFast();
3107 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3109 for (Int_t i=0; i<entries; i++){
3110 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3111 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3112 if (!kalmanP) continue;
3113 if (!kalmanM) continue;
3115 AliRelAlignerKalman *kalman = NULL;
3116 if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3117 kalman = new AliRelAlignerKalman(*kalmanP);
3118 kalman->Merge(kalmanM);
3120 else if (kalmanP->GetRunNumber() == 0) {
3121 kalman = new AliRelAlignerKalman(*kalmanM);
3123 else if (kalmanM->GetRunNumber() == 0) {
3124 kalman = new AliRelAlignerKalman(*kalmanP);
3129 array->AddAt(kalman,i);