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 "AliTPCcalibDB.h"
44 #include "AliTPCCalPad.h"
45 #include "AliTPCCalROC.h"
46 #include "AliTPCROC.h"
47 #include "AliTPCmapper.h"
48 #include "AliTPCParam.h"
49 #include "AliTPCCalibRaw.h"
50 #include "AliTPCPreprocessorOnline.h"
51 #include "AliTPCdataQA.h"
53 #include "AliTPCcalibDButil.h"
54 #include "AliTPCCalibVdrift.h"
55 #include "AliMathBase.h"
57 ClassImp(AliTPCcalibDButil)
58 AliTPCcalibDButil::AliTPCcalibDButil() :
66 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
78 fRefPedestalMasked(0x0),
82 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
83 fRefPulserMasked(0x0),
90 fRefALTROAcqStart(0x0),
91 fRefALTROAcqStop(0x0),
96 fMapper(new AliTPCmapper(0x0)),
100 fPulTmaxLimitAbs(1.5),
103 fRuns(0), // run list with OCDB info
104 fRunsStart(0), // start time for given run
105 fRunsStop(0) // stop time for given run
111 //_____________________________________________________________________________________
112 AliTPCcalibDButil::~AliTPCcalibDButil()
117 delete fPulserOutlier;
118 delete fRefPulserOutlier;
120 if (fRefPadNoise) delete fRefPadNoise;
121 if (fRefPedestals) delete fRefPedestals;
122 if (fRefPedestalMasked) delete fRefPedestalMasked;
123 if (fRefPulserTmean) delete fRefPulserTmean;
124 if (fRefPulserTrms) delete fRefPulserTrms;
125 if (fRefPulserQmean) delete fRefPulserQmean;
126 if (fRefPulserMasked) delete fRefPulserMasked;
127 if (fRefCETmean) delete fRefCETmean;
128 if (fRefCETrms) delete fRefCETrms;
129 if (fRefCEQmean) delete fRefCEQmean;
130 if (fRefCEMasked) delete fRefCEMasked;
131 if (fRefALTROFPED) delete fRefALTROFPED;
132 if (fRefALTROZsThr) delete fRefALTROZsThr;
133 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
134 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
135 if (fRefALTROMasked) delete fRefALTROMasked;
136 if (fRefCalibRaw) delete fRefCalibRaw;
137 if (fCurrentRefMap) delete fCurrentRefMap;
139 //_____________________________________________________________________________________
140 void AliTPCcalibDButil::UpdateFromCalibDB()
143 // Update pointers from calibDB
145 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
146 fPadNoise=fCalibDB->GetPadNoise();
147 fPedestals=fCalibDB->GetPedestals();
148 fPulserTmean=fCalibDB->GetPulserTmean();
149 fPulserTrms=fCalibDB->GetPulserTrms();
150 fPulserQmean=fCalibDB->GetPulserQmean();
151 fCETmean=fCalibDB->GetCETmean();
152 fCETrms=fCalibDB->GetCETrms();
153 fCEQmean=fCalibDB->GetCEQmean();
154 fALTROMasked=fCalibDB->GetALTROMasked();
155 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
156 fCalibRaw=fCalibDB->GetCalibRaw();
157 fDataQA=fCalibDB->GetDataQA();
158 UpdatePulserOutlierMap();
160 UpdateRefDataFromOCDB();
162 //_____________________________________________________________________________________
163 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
164 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
167 // Process the CE data for this run
168 // the return TVectorD arrays contian the results of the fit
169 // noutliersCE contains the number of pads marked as outliers,
170 // not including masked and edge pads
173 //retrieve CE and ALTRO data
175 TString fitString(fitFormula);
176 fitString.ReplaceAll("++","#");
177 Int_t ndim=fitString.CountChar('#')+2;
178 fitResultsA.ResizeTo(ndim);
179 fitResultsC.ResizeTo(ndim);
188 if (outCE) out=outCE;
189 else out=new AliTPCCalPad("outCE","outCE");
190 AliTPCCalROC *rocMasked=0x0;
191 //loop over all channels
192 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
193 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
194 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
195 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
197 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
201 //add time offset to IROCs
202 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
203 rocData->Add(fIrocTimeOffset);
205 UInt_t nrows=rocData->GetNrows();
206 for (UInt_t irow=0;irow<nrows;++irow){
207 UInt_t npads=rocData->GetNPads(irow);
208 for (UInt_t ipad=0;ipad<npads;++ipad){
209 rocOut->SetValue(irow,ipad,0);
210 //exclude masked pads
211 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
212 rocOut->SetValue(irow,ipad,1);
215 //exclude first two rows in IROC and last two rows in OROC
217 if (irow<2) rocOut->SetValue(irow,ipad,1);
219 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
222 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
223 Float_t valTmean=rocData->GetValue(irow,ipad);
224 //exclude values that are exactly 0
226 rocOut->SetValue(irow,ipad,1);
229 // exclude channels with too large variations
230 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
231 rocOut->SetValue(irow,ipad,1);
239 Float_t chi2Af,chi2Cf;
240 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
243 if (!outCE) delete out;
245 //_____________________________________________________________________________________
246 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
247 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
248 Float_t &driftTimeA, Float_t &driftTimeC )
251 // Calculate statistical information from the CE graphs for drift time and charge
255 vecTEntries.ResizeTo(72);
256 vecTMean.ResizeTo(72);
257 vecTRMS.ResizeTo(72);
258 vecTMedian.ResizeTo(72);
259 vecQEntries.ResizeTo(72);
260 vecQMean.ResizeTo(72);
261 vecQRMS.ResizeTo(72);
262 vecQMedian.ResizeTo(72);
273 TObjArray *arrT=fCalibDB->GetCErocTtime();
274 TObjArray *arrQ=fCalibDB->GetCErocQtime();
276 for (Int_t isec=0;isec<74;++isec){
277 TGraph *gr=(TGraph*)arrT->At(isec);
280 Int_t npoints = gr->GetN();
281 values.ResizeTo(npoints);
283 //skip first points, theres always some problems with finding the CE position
284 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
285 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
286 values[nused]=gr->GetY()[ipoint];
291 if (isec<72) vecTEntries[isec]= nused;
294 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
295 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
296 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
297 } else if (isec==72){
298 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
299 } else if (isec==73){
300 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
306 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
307 TGraph *gr=(TGraph*)arrQ->At(isec);
310 Int_t npoints = gr->GetN();
311 values.ResizeTo(npoints);
313 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
314 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
315 values[nused]=gr->GetY()[ipoint];
320 vecQEntries[isec]= nused;
322 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
323 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
324 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
330 //_____________________________________________________________________________________
331 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
332 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
333 Int_t &nonMaskedZero)
336 // process noise data
337 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
338 // OROCs small pads [2] and OROCs large pads [3]
339 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
340 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
343 //set proper size and reset
344 const UInt_t infoSize=4;
345 vNoiseMean.ResizeTo(infoSize);
346 vNoiseMeanSenRegions.ResizeTo(infoSize);
347 vNoiseRMS.ResizeTo(infoSize);
348 vNoiseRMSSenRegions.ResizeTo(infoSize);
350 vNoiseMeanSenRegions.Zero();
352 vNoiseRMSSenRegions.Zero();
355 TVectorD c(infoSize);
356 TVectorD cs(infoSize);
360 //retrieve noise and ALTRO data
361 if (!fPadNoise) return;
362 AliTPCCalROC *rocMasked=0x0;
363 //create IROC, OROC1, OROC2 and sensitive region masks
364 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
365 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
366 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
367 UInt_t nrows=noiseROC->GetNrows();
368 for (UInt_t irow=0;irow<nrows;++irow){
369 UInt_t npads=noiseROC->GetNPads(irow);
370 for (UInt_t ipad=0;ipad<npads;++ipad){
371 //don't use masked channels;
372 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
373 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
380 if ( !(noiseVal<10000000) ){
381 printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
384 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
385 Int_t masksen=1; // sensitive pards are not masked (0)
386 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
387 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
389 if (irow>19&&irow<46){
390 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
393 vNoiseMean[type]+=noiseVal;
394 vNoiseRMS[type]+=noiseVal*noiseVal;
397 vNoiseMeanSenRegions[type]+=noiseVal;
398 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
403 //define sensive regions
404 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
406 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
407 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
409 if ((Int_t)irow<par.GetNRowUp1()){
412 vNoiseMean[type]+=noiseVal;
413 vNoiseRMS[type]+=noiseVal*noiseVal;
416 vNoiseMeanSenRegions[type]+=noiseVal;
417 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
423 vNoiseMean[type]+=noiseVal;
424 vNoiseRMS[type]+=noiseVal*noiseVal;
427 vNoiseMeanSenRegions[type]+=noiseVal;
428 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
435 vNoiseMean[type]+=noiseVal;
436 vNoiseRMS[type]+=noiseVal*noiseVal;
439 vNoiseMeanSenRegions[type]+=noiseVal;
440 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
445 }//end loop sectors (rocs)
447 //calculate mean and RMS
448 const Double_t verySmall=0.0000000001;
449 for (UInt_t i=0;i<infoSize;++i){
456 // printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
457 mean=vNoiseMean[i]/c[i];
459 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
464 if (cs[i]>verySmall){
465 meanSen=vNoiseMeanSenRegions[i]/cs[i];
466 rmsSen=vNoiseRMSSenRegions[i];
467 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
469 vNoiseMeanSenRegions[i]=meanSen;
470 vNoiseRMSSenRegions[i]=rmsSen;
474 //_____________________________________________________________________________________
475 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
478 // Process the Pulser information
479 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
482 const UInt_t infoSize=4;
483 //reset counters to error number
484 vMeanTime.ResizeTo(infoSize);
487 TVectorD c(infoSize);
488 //retrieve pulser and ALTRO data
489 if (!fPulserTmean) return;
492 AliTPCCalROC *rocOut=0x0;
493 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
494 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
495 if (!tmeanROC) continue;
496 rocOut=fPulserOutlier->GetCalROC(isec);
497 UInt_t nchannels=tmeanROC->GetNchannels();
498 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
499 if (rocOut && rocOut->GetValue(ichannel)) continue;
500 Float_t val=tmeanROC->GetValue(ichannel);
502 vMeanTime[type]+=val;
507 for (UInt_t itype=0; itype<infoSize; ++itype){
508 if (c[itype]>0) vMeanTime[itype]/=c[itype];
509 else vMeanTime[itype]=0;
512 //_____________________________________________________________________________________
513 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
516 // Get Values from ALTRO configuration data
519 if (!fALTROMasked) return;
521 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
522 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
523 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
524 if (rocMasked->GetValue(ichannel)) ++nMasked;
528 //_____________________________________________________________________________________
529 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
532 // Proces Goofie values, return statistical information of the currently set goofieArray
533 // The meaning of the entries are given below
535 1 TPC_ANODE_I_A00_STAT
537 3 TPC_DVM_DriftVelocity
542 8 TPC_DVM_NumberOfSparks
543 9 TPC_DVM_PeakAreaFar
544 10 TPC_DVM_PeakAreaNear
545 11 TPC_DVM_PeakPosFar
546 12 TPC_DVM_PeakPosNear
552 18 TPC_DVM_TemperatureS1
556 vecEntries.ResizeTo(nsensors);
557 vecMedian.ResizeTo(nsensors);
558 vecMean.ResizeTo(nsensors);
559 vecRMS.ResizeTo(nsensors);
566 Double_t kEpsilon=0.0000000001;
567 Double_t kBig=100000000000.;
568 Int_t nsensors = fGoofieArray->NumSensors();
569 vecEntries.ResizeTo(nsensors);
570 vecMedian.ResizeTo(nsensors);
571 vecMean.ResizeTo(nsensors);
572 vecRMS.ResizeTo(nsensors);
574 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
575 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
576 if (gsensor && gsensor->GetGraph()){
577 Int_t npoints = gsensor->GetGraph()->GetN();
579 values.ResizeTo(npoints);
581 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
582 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
583 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
584 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
589 vecEntries[isensor]= nused;
591 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
592 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
593 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
598 //_____________________________________________________________________________________
599 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
602 // check the variations of the pedestal data to the reference pedestal data
603 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
606 TVectorF vThres(npar); //thresholds
607 Int_t nActive=0; //number of active channels
609 //reset and set thresholds
610 pedestalDeviations.ResizeTo(npar);
611 for (Int_t i=0;i<npar;++i){
612 pedestalDeviations.GetMatrixArray()[i]=0;
613 vThres.GetMatrixArray()[i]=(i+1)*.5;
615 //check all needed data is available
616 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
617 //loop over all channels
618 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
619 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
620 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
621 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
622 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
623 UInt_t nrows=mROC->GetNrows();
624 for (UInt_t irow=0;irow<nrows;++irow){
625 UInt_t npads=mROC->GetNPads(irow);
626 for (UInt_t ipad=0;ipad<npads;++ipad){
627 //don't use masked channels;
628 if (mROC ->GetValue(irow,ipad)) continue;
629 if (mRefROC->GetValue(irow,ipad)) continue;
630 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
631 for (Int_t i=0;i<npar;++i){
632 if (deviation>vThres[i])
633 ++pedestalDeviations.GetMatrixArray()[i];
640 for (Int_t i=0;i<npar;++i){
641 pedestalDeviations.GetMatrixArray()[i]/=nActive;
645 //_____________________________________________________________________________________
646 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
649 // check the variations of the noise data to the reference noise data
650 // thresholds are 5, 10, 15 and 20 percent respectively.
653 TVectorF vThres(npar); //thresholds
654 Int_t nActive=0; //number of active channels
656 //reset and set thresholds
657 noiseDeviations.ResizeTo(npar);
658 for (Int_t i=0;i<npar;++i){
659 noiseDeviations.GetMatrixArray()[i]=0;
660 vThres.GetMatrixArray()[i]=(i+1)*.05;
662 //check all needed data is available
663 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
664 //loop over all channels
665 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
666 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
667 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
668 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
669 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
670 UInt_t nrows=mROC->GetNrows();
671 for (UInt_t irow=0;irow<nrows;++irow){
672 UInt_t npads=mROC->GetNPads(irow);
673 for (UInt_t ipad=0;ipad<npads;++ipad){
674 //don't use masked channels;
675 if (mROC ->GetValue(irow,ipad)) continue;
676 if (mRefROC->GetValue(irow,ipad)) continue;
677 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
678 for (Int_t i=0;i<npar;++i){
679 if (deviation>vThres[i])
680 ++noiseDeviations.GetMatrixArray()[i];
687 for (Int_t i=0;i<npar;++i){
688 noiseDeviations.GetMatrixArray()[i]/=nActive;
692 //_____________________________________________________________________________________
693 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
694 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
697 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
698 // thresholds are .5, 1, 5 and 10 percent respectively.
702 TVectorF vThres(npar); //thresholds
703 Int_t nActive=0; //number of active channels
705 //reset and set thresholds
706 pulserQdeviations.ResizeTo(npar);
707 for (Int_t i=0;i<npar;++i){
708 pulserQdeviations.GetMatrixArray()[i]=0;
713 vThres.GetMatrixArray()[0]=.005;
714 vThres.GetMatrixArray()[1]=.01;
715 vThres.GetMatrixArray()[2]=.05;
716 vThres.GetMatrixArray()[3]=.1;
717 //check all needed data is available
718 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
720 UpdateRefPulserOutlierMap();
721 //loop over all channels
722 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
723 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
724 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
725 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
726 // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
727 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
728 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
729 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
730 Float_t pt_mean=ptROC->GetMean(oROC);
731 UInt_t nrows=mROC->GetNrows();
732 for (UInt_t irow=0;irow<nrows;++irow){
733 UInt_t npads=mROC->GetNPads(irow);
734 for (UInt_t ipad=0;ipad<npads;++ipad){
735 //don't use masked channels;
736 if (mROC ->GetValue(irow,ipad)) continue;
737 if (mRefROC->GetValue(irow,ipad)) continue;
738 //don't user edge pads
739 if (ipad==0||ipad==npads-1) continue;
741 Float_t pq=pqROC->GetValue(irow,ipad);
742 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
743 Float_t pt=ptROC->GetValue(irow,ipad);
744 // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
746 Float_t deviation=TMath::Abs(pq/pqRef-1);
747 for (Int_t i=0;i<npar;++i){
748 if (deviation>vThres[i])
749 ++pulserQdeviations.GetMatrixArray()[i];
751 if (pqRef>11&&pq<11) ++npadsOffAdd;
754 if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
760 for (Int_t i=0;i<npar;++i){
761 pulserQdeviations.GetMatrixArray()[i]/=nActive;
766 //_____________________________________________________________________________________
767 void AliTPCcalibDButil::UpdatePulserOutlierMap()
772 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
774 //_____________________________________________________________________________________
775 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
780 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
782 //_____________________________________________________________________________________
783 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
786 // Create a map that contains outliers from the Pulser calibration data.
787 // The outliers include masked channels, edge pads and pads with
788 // too large timing and charge variations.
789 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
790 // those do not contain masked and edge pads
794 pulOut->Multiply(0.);
798 AliTPCCalROC *rocMasked=0x0;
802 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
803 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
804 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
805 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
806 if (!tmeanROC||!qmeanROC) {
807 //reset outliers in this ROC
808 outROC->Multiply(0.);
811 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
813 // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
814 // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
815 UInt_t nrows=tmeanROC->GetNrows();
816 for (UInt_t irow=0;irow<nrows;++irow){
817 UInt_t npads=tmeanROC->GetNPads(irow);
818 for (UInt_t ipad=0;ipad<npads;++ipad){
819 Int_t outlier=0,masked=0;
820 Float_t q=qmeanROC->GetValue(irow,ipad);
821 Float_t t=tmeanROC->GetValue(irow,ipad);
822 //masked channels are outliers
823 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
824 //edge pads are outliers
825 if (ipad==0||ipad==npads-1) masked=1;
826 //channels with too large charge or timing deviation from the meadian are outliers
827 // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
828 if (q<fPulQminLimit && !masked) outlier=1;
830 if ( !(q<10000000) || !(t<10000000)) outlier=1;
831 outROC->SetValue(irow,ipad,outlier+masked);
832 fNpulserOutliers+=outlier;
837 //_____________________________________________________________________________________
838 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
841 // Create pad time0 object from pulser and/or CE data, depending on the selected model
842 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
843 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
844 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
846 // In case model 2 is invoked - gy arival time gradient is also returned
850 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
851 // decide between different models
852 if (model==0||model==1){
854 if (model==1) ProcessPulser(vMean);
855 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
856 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
857 if (!rocPulTmean) continue;
858 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
859 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
860 Float_t mean=rocPulTmean->GetMean(rocOut);
861 //treat case where a whole partition is masked
862 if (mean==0) mean=rocPulTmean->GetMean();
867 UInt_t nrows=rocTime0->GetNrows();
868 for (UInt_t irow=0;irow<nrows;++irow){
869 UInt_t npads=rocTime0->GetNPads(irow);
870 for (UInt_t ipad=0;ipad<npads;++ipad){
871 Float_t time=rocPulTmean->GetValue(irow,ipad);
872 //in case of an outlier pad use the mean of the altro values.
873 //This should be the most precise guess in that case.
874 if (rocOut->GetValue(irow,ipad)) {
875 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
876 if (time==0) time=mean;
878 Float_t val=time-mean;
879 rocTime0->SetValue(irow,ipad,val);
883 } else if (model==2){
884 Double_t pgya,pgyc,pchi2a,pchi2c;
885 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
886 fCETmean->Add(padPulser,-1.);
888 AliTPCCalPad outCE("outCE","outCE");
890 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
891 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
892 // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
893 if (!padFit) { delete padPulser; return 0;}
896 fCETmean->Add(padPulser,1.);
897 padTime0->Add(fCETmean);
898 padTime0->Add(padFit,-1);
903 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
904 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
905 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
906 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
907 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
908 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
909 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
910 Float_t mean=rocPulTmean->GetMean(rocOutPul);
911 if (mean==0) mean=rocPulTmean->GetMean();
912 UInt_t nrows=rocTime0->GetNrows();
913 for (UInt_t irow=0;irow<nrows;++irow){
914 UInt_t npads=rocTime0->GetNPads(irow);
915 for (UInt_t ipad=0;ipad<npads;++ipad){
916 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
917 if (rocOutCE->GetValue(irow,ipad)){
918 Float_t valOut=rocCEfit->GetValue(irow,ipad);
919 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
920 rocTime0->SetValue(irow,ipad,valOut);
928 Double_t median = padTime0->GetMedian();
929 padTime0->Add(-median); // normalize to median
932 //_____________________________________________________________________________________
933 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
935 if (roc==0) return 0.;
936 const Int_t sector=roc->GetSector();
937 AliTPCROC *tpcRoc=AliTPCROC::Instance();
938 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
942 //loop over a small range around the requested pad (+-10 rows/pads)
943 for (Int_t irow=row-10;irow<row+10;++irow){
944 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
945 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
946 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
947 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
948 if (altroRoc!=altroCurr) continue;
949 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
950 Float_t val=roc->GetValue(irow,ipad);
958 //_____________________________________________________________________________________
959 void AliTPCcalibDButil::SetRefFile(const char* filename)
962 // load cal pad objects form the reference file
964 TDirectory *currDir=gDirectory;
966 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
967 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
969 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
970 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
971 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
973 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
974 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
975 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
977 // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
978 // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
979 // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
980 // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
981 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
985 //_____________________________________________________________________________________
986 void AliTPCcalibDButil::UpdateRefDataFromOCDB()
989 // set reference data from OCDB Reference map
992 AliWarning("Referenc map not set!");
997 AliCDBEntry* entry = 0x0;
998 Bool_t hasAnyChanged=kFALSE;
1001 cdbPath="TPC/Calib/Pedestals";
1002 if (HasRefChanged(cdbPath.Data())){
1003 hasAnyChanged=kTRUE;
1004 //delete old entries
1005 if (fRefPedestals) delete fRefPedestals;
1006 if (fRefPedestalMasked) delete fRefPedestalMasked;
1007 fRefPedestals=fRefPedestalMasked=0x0;
1009 entry=GetRefEntry(cdbPath.Data());
1011 entry->SetOwner(kTRUE);
1012 fRefPedestals=GetRefCalPad(entry);
1014 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1019 cdbPath="TPC/Calib/PadNoise";
1020 if (HasRefChanged(cdbPath.Data())){
1021 hasAnyChanged=kTRUE;
1023 if (fRefPadNoise) delete fRefPadNoise;
1026 entry=GetRefEntry(cdbPath.Data());
1028 entry->SetOwner(kTRUE);
1029 fRefPadNoise=GetRefCalPad(entry);
1035 cdbPath="TPC/Calib/Pulser";
1036 if (HasRefChanged(cdbPath.Data())){
1037 hasAnyChanged=kTRUE;
1038 //delete old entries
1039 if (fRefPulserTmean) delete fRefPulserTmean;
1040 if (fRefPulserTrms) delete fRefPulserTrms;
1041 if (fRefPulserQmean) delete fRefPulserQmean;
1042 if (fRefPulserMasked) delete fRefPulserMasked;
1043 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1045 entry=GetRefEntry(cdbPath.Data());
1047 entry->SetOwner(kTRUE);
1048 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1049 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1050 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1052 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1057 cdbPath="TPC/Calib/CE";
1058 if (HasRefChanged(cdbPath.Data())){
1059 hasAnyChanged=kTRUE;
1060 //delete old entries
1061 if (fRefCETmean) delete fRefCETmean;
1062 if (fRefCETrms) delete fRefCETrms;
1063 if (fRefCEQmean) delete fRefCEQmean;
1064 if (fRefCEMasked) delete fRefCEMasked;
1065 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1067 entry=GetRefEntry(cdbPath.Data());
1069 entry->SetOwner(kTRUE);
1070 fRefCETmean=GetRefCalPad(entry,"CETmean");
1071 fRefCETrms=GetRefCalPad(entry,"CETrms");
1072 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1074 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1079 cdbPath="TPC/Calib/AltroConfig";
1080 if (HasRefChanged(cdbPath.Data())){
1081 hasAnyChanged=kTRUE;
1082 //delete old entries
1083 if (fRefALTROFPED) delete fRefALTROFPED;
1084 if (fRefALTROZsThr) delete fRefALTROZsThr;
1085 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1086 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1087 if (fRefALTROMasked) delete fRefALTROMasked;
1088 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1090 entry=GetRefEntry(cdbPath.Data());
1092 entry->SetOwner(kTRUE);
1093 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1094 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1095 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1096 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1097 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1104 cdbPath="TPC/Calib/Raw";
1105 if (HasRefChanged(cdbPath.Data())){
1106 hasAnyChanged=kTRUE;
1108 if (fRefCalibRaw) delete fRefCalibRaw;
1110 entry=GetRefEntry(cdbPath.Data());
1112 entry->SetOwner(kTRUE);
1113 TObjArray *arr=(TObjArray*)entry->GetObject();
1115 AliError(Form("Could not get object from entry '%s'\nPlese check!!!",entry->GetId().GetPath().Data()));
1117 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1124 cdbPath="TPC/Calib/QA";
1125 if (HasRefChanged(cdbPath.Data())){
1126 hasAnyChanged=kTRUE;
1128 if (fRefDataQA) delete fRefDataQA;
1130 entry=GetRefEntry(cdbPath.Data());
1132 entry->SetOwner(kTRUE);
1133 fDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1135 AliError(Form("Could not get object from entry '%s'\nPlese check!!!",entry->GetId().GetPath().Data()));
1137 fRefDataQA=(AliTPCdataQA*)fDataQA->Clone();
1144 //update current reference maps
1146 if (fCurrentRefMap) delete fCurrentRefMap;
1147 fCurrentRefMap=(TMap*)fRefMap->Clone();
1150 //_____________________________________________________________________________________
1151 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1154 // TObjArray object type case
1155 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1157 AliTPCCalPad *pad=0x0;
1158 TObjArray *arr=(TObjArray*)entry->GetObject();
1160 AliError(Form("Could not get object from entry '%s'\nPlese check!!!",entry->GetId().GetPath().Data()));
1163 pad=(AliTPCCalPad*)arr->FindObject(objName);
1165 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlese check!!!",objName,entry->GetId().GetPath().Data()));
1168 return (AliTPCCalPad*)pad->Clone();
1170 //_____________________________________________________________________________________
1171 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1174 // AliTPCCalPad object type case
1175 // cast object to a calPad and store it in 'pad'
1177 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1179 AliError(Form("Could not get object from entry '%s'\nPlese check!!!",entry->GetId().GetPath().Data()));
1182 pad=(AliTPCCalPad*)pad->Clone();
1185 //_____________________________________________________________________________________
1186 AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1189 // set altro masked channel map for 'cdbPath'
1191 AliTPCCalPad* pad=0x0;
1192 const Int_t run=GetReferenceRun(cdbPath);
1194 AliError(Form("Could not get reference run number for object '%s'\nPlese check availability!!!",cdbPath));
1197 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1199 AliError(Form("Could not get reference object '%s'\nPlese check availability!!!",cdbPath));
1202 pad=GetRefCalPad(entry,"Masked");
1203 if (pad) pad->SetNameTitle(name,name);
1204 entry->SetOwner(kTRUE);
1208 //_____________________________________________________________________________________
1209 void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1211 // Get Reference map
1213 if (run<0) run=fCalibDB->GetRun();
1214 TString cdbPath="TPC/Calib/Ref";
1215 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1217 AliError(Form("Could not get reference object '%s'\nPlese check availability!!!",cdbPath.Data()));
1221 entry->SetOwner(kTRUE);
1222 fRefMap=(TMap*)(entry->GetObject());
1223 AliCDBId &id=entry->GetId();
1224 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1226 //_____________________________________________________________________________________
1227 Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1230 // check whether a reference cdb entry has changed
1232 if (!fCurrentRefMap) return kTRUE;
1233 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1236 //_____________________________________________________________________________________
1237 AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1240 // get the reference AliCDBEntry for 'cdbPath'
1242 const Int_t run=GetReferenceRun(cdbPath);
1244 AliError(Form("Could not get reference run number for object '%s'\nPlese check availability!!!",cdbPath));
1247 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1249 AliError(Form("Could not get reference object '%s'\nPlese check availability!!!",cdbPath));
1254 //_____________________________________________________________________________________
1255 const Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type){
1257 // Get reference run number for the specified OCDB path
1259 if (!fCurrentRefMap) return -2;
1260 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1261 if (!str) return -2;
1262 return str->GetString().Atoi();
1264 //_____________________________________________________________________________________
1265 const Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1267 // Get reference run number for the specified OCDB path
1269 if (!fRefMap) return -1;
1270 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1271 if (!str) return -1;
1272 return str->GetString().Atoi();
1274 //_____________________________________________________________________________________
1275 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
1277 // Author: marian.ivanov@cern.ch
1279 // Create outlier map for CE study
1281 // Return value - outlyer map
1282 // noutlyersCE - number of outlyers
1283 // minSignal - minimal total Q signal
1284 // cutRMSMin - minimal width of the signal in respect to the median
1285 // cutRMSMax - maximal width of the signal in respect to the median
1286 // cutMaxDistT - maximal deviation from time median per chamber
1288 // Outlyers criteria:
1289 // 0. Exclude masked pads
1290 // 1. Exclude first two rows in IROC and last two rows in OROC
1291 // 2. Exclude edge pads
1292 // 3. Exclude channels with too large variations
1293 // 4. Exclude pads with too small signal
1294 // 5. Exclude signal with outlyers RMS
1295 // 6. Exclude channels to far from the chamber median
1297 //create outlier map
1298 AliTPCCalPad *out=ceOut;
1299 if (!out) out= new AliTPCCalPad("outCE","outCE");
1300 AliTPCCalROC *rocMasked=0x0;
1301 if (!fCETmean) return 0;
1302 if (!fCETrms) return 0;
1303 if (!fCEQmean) return 0;
1305 //loop over all channels
1307 Double_t rmsMedian = fCETrms->GetMedian();
1308 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1309 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1310 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1311 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1312 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1313 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1314 Double_t trocMedian = rocData->GetMedian();
1317 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1323 UInt_t nrows=rocData->GetNrows();
1324 for (UInt_t irow=0;irow<nrows;++irow){
1325 UInt_t npads=rocData->GetNPads(irow);
1326 for (UInt_t ipad=0;ipad<npads;++ipad){
1327 rocOut->SetValue(irow,ipad,0);
1328 Float_t valTmean=rocData->GetValue(irow,ipad);
1329 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1330 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1331 //0. exclude masked pads
1332 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1333 rocOut->SetValue(irow,ipad,1);
1336 //1. exclude first two rows in IROC and last two rows in OROC
1338 if (irow<2) rocOut->SetValue(irow,ipad,1);
1340 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1342 //2. exclude edge pads
1343 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1344 //exclude values that are exactly 0
1346 rocOut->SetValue(irow,ipad,1);
1349 //3. exclude channels with too large variations
1350 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1351 rocOut->SetValue(irow,ipad,1);
1355 //4. exclude channels with too small signal
1356 if (valQmean<minSignal) {
1357 rocOut->SetValue(irow,ipad,1);
1361 //5. exclude channels with too small rms
1362 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1363 rocOut->SetValue(irow,ipad,1);
1367 //6. exclude channels to far from the chamber median
1368 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1369 rocOut->SetValue(irow,ipad,1);
1380 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1382 // Author: marian.ivanov@cern.ch
1384 // Create outlier map for Pulser
1386 // Return value - outlyer map
1387 // noutlyersPulser - number of outlyers
1388 // cutTime - absolute cut - distance to the median of chamber
1389 // cutnRMSQ - nsigma cut from median q distribution per chamber
1390 // cutnRMSrms - nsigma cut from median rms distribution
1391 // Outlyers criteria:
1392 // 0. Exclude masked pads
1393 // 1. Exclude time outlyers (default 3 time bins)
1394 // 2. Exclude q outlyers (default 5 sigma)
1395 // 3. Exclude rms outlyers (default 5 sigma)
1397 AliTPCCalPad *out=pulserOut;
1398 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1399 AliTPCCalROC *rocMasked=0x0;
1400 if (!fPulserTmean) return 0;
1401 if (!fPulserTrms) return 0;
1402 if (!fPulserQmean) return 0;
1404 //loop over all channels
1406 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1407 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1408 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1409 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1410 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1411 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1413 Double_t rocMedianT = rocData->GetMedian();
1414 Double_t rocMedianQ = rocPulserQ->GetMedian();
1415 Double_t rocRMSQ = rocPulserQ->GetRMS();
1416 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1417 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1418 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1419 rocOut->SetValue(ichannel,0);
1420 Float_t valTmean=rocData->GetValue(ichannel);
1421 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1422 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1424 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1425 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1426 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1427 rocOut->SetValue(ichannel,isOut);
1428 if (isOut) noutliersPulser++;
1435 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1437 // Author : Marian Ivanov
1438 // Create pad time0 correction map using information from the CE and from pulser
1441 // Return PadTime0 to be used for time0 relative alignment
1442 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1443 // using $ALICE_ROOT/TPC/script/gui application
1445 // fitResultsA - fitParameters A side
1446 // fitResultsC - fitParameters C side
1447 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1448 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1452 // 1. Find outlier map for CE
1453 // 2. Find outlier map for Pulser
1454 // 3. Replace outlier by median at given sector (median without outliers)
1455 // 4. Substract from the CE data pulser
1456 // 5. Fit the CE with formula
1457 // 5.1) (IROC-OROC) offset
1461 // 5.5) (IROC-OROC)*(lx-xmid)
1463 // 6. Substract gy fit dependence from the CE data
1464 // 7. Add pulser back to CE data
1465 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1466 // 9. return CE data
1468 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1469 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1472 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)";
1473 // output for fit formula
1474 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)";
1475 // gy part of formula
1476 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)";
1479 if (!fCETmean) return 0;
1480 Double_t pgya,pgyc,pchi2a,pchi2c;
1481 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1482 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1484 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1485 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1486 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1487 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1488 padPulser->SetName("padPulser");
1489 padPulserOut->SetName("padPulserOut");
1490 padCE->SetName("padCE");
1491 padCEIn->SetName("padCEIn");
1492 padCEOut->SetName("padCEOut");
1493 padOut->SetName("padOut");
1496 // make combined outlyers map
1497 // and replace outlyers in maps with median for chamber
1499 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1500 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1501 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1502 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1503 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1504 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1505 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1506 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1507 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1508 if (rocPulserOut->GetValue(ichannel)>0) {
1509 rocPulser->SetValue(ichannel,pulserMedian);
1510 rocOut->SetValue(ichannel,1);
1512 if (rocCEOut->GetValue(ichannel)>0) {
1513 rocCE->SetValue(ichannel,ceMedian);
1514 rocOut->SetValue(ichannel,1);
1519 // remove pulser time 0
1521 padCE->Add(padPulser,-1);
1526 Float_t chi2Af,chi2Cf;
1527 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1531 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1532 padCEFitGY->SetName("padCEFitGy");
1534 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1535 padCEFit->SetName("padCEFit");
1537 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1538 padCEDiff->SetName("padCEDiff");
1539 padCEDiff->Add(padCEFit,-1.);
1542 padCE->Add(padCEFitGY,-1.);
1544 padCE->Add(padPulser,1.);
1545 Double_t padmedian = padCE->GetMedian();
1546 padCE->Add(-padmedian); // normalize to median
1548 // Replace outliers by fit value - median of diff per given chamber -GY fit
1550 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1551 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1552 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1553 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1554 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1555 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1557 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1558 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1559 if (rocOut->GetValue(ichannel)==0) continue;
1560 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1561 rocCE->SetValue(ichannel,value);
1567 //dump to the file - result can be visualized
1568 AliTPCPreprocessorOnline preprocesor;
1569 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1570 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1571 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1572 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1574 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1575 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1577 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1578 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1579 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1580 preprocesor.DumpToFile(dumpfile);
1583 delete padPulserOut;
1596 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1598 // find the closest point to xref in x direction
1599 // return dx and value
1601 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1602 if (index<0) index=0;
1603 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1604 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1605 dx = xref-graph->GetX()[index];
1606 y = graph->GetY()[index];
1611 Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1613 // Get the correction of the trigger offset
1614 // combining information from the laser track calibration
1615 // and from cosmic calibration
1618 // timeStamp - tim stamp in seconds
1619 // deltaT - integration period to calculate offset
1620 // deltaTLaser -max validity of laser data
1621 // valType - 0 - median, 1- mean
1623 // Integration vaues are just recomendation - if not possible to get points
1624 // automatically increase the validity by factor 2
1625 // (recursive algorithm until one month of data taking)
1628 const Float_t kLaserCut=0.0005;
1629 const Int_t kMaxPeriod=3600*24*30*3; // 3 month max
1630 const Int_t kMinPoints=20;
1632 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1634 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1636 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1637 if (!array) return 0;
1639 TGraphErrors *laserA[3]={0,0,0};
1640 TGraphErrors *laserC[3]={0,0,0};
1641 TGraphErrors *cosmicAll=0;
1642 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1643 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1644 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1647 if (!cosmicAll) return 0;
1648 Int_t nmeasC=cosmicAll->GetN();
1649 Float_t *tdelta = new Float_t[nmeasC];
1651 for (Int_t i=0;i<nmeasC;i++){
1652 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1653 Float_t ccosmic=cosmicAll->GetY()[i];
1654 Double_t yA=0,yC=0,dA=0,dC=0;
1655 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1656 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1657 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1658 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1660 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1662 if (TMath::Abs(yA-yC)<kLaserCut) {
1665 if (i%2==0) claser=yA;
1666 if (i%2==1) claser=yC;
1668 tdelta[nused]=ccosmic-claser;
1671 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1672 Double_t median = TMath::Median(nused,tdelta);
1673 Double_t mean = TMath::Mean(nused,tdelta);
1675 return (valType==0) ? median:mean;
1678 Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1680 // Get the correction of the drift velocity
1681 // combining information from the laser track calibration
1682 // and from cosmic calibration
1684 // dist - return value - distance to closest point in graph
1686 // timeStamp - tim stamp in seconds
1687 // deltaT - integration period to calculate time0 offset
1688 // deltaTLaser -max validity of laser data
1689 // valType - 0 - median, 1- mean
1691 // Integration vaues are just recomendation - if not possible to get points
1692 // automatically increase the validity by factor 2
1693 // (recursive algorithm until one month of data taking)
1697 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1699 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1701 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1702 if (!array) return 0;
1703 TGraphErrors *cosmicAll=0;
1704 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1705 if (!cosmicAll) return 0;
1707 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1709 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1710 Double_t vcosmic= AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1711 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1712 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1719 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1720 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1721 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1723 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1724 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1725 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1726 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1728 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1729 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1735 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1738 // Create a default name for the gui file
1741 return Form("guiRefTreeRun%s.root",GetRefValidity());
1744 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1747 // Create a gui reference tree
1748 // if dirname and filename are empty default values will be used
1749 // this is the recommended way of using this function
1750 // it allows to check whether a file with the given run validity alredy exists
1752 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1753 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1757 TString file=filename;
1758 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1760 AliTPCPreprocessorOnline prep;
1761 //noise and pedestals
1762 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1763 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1764 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1766 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1767 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1768 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1769 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1771 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1772 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1773 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1774 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1776 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1777 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1778 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1779 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1780 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1782 AliTPCdataQA *dataQA=fRefDataQA;
1784 if (dataQA->GetNLocalMaxima())
1785 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1786 if (dataQA->GetMaxCharge())
1787 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1788 if (dataQA->GetMeanCharge())
1789 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1790 if (dataQA->GetNoThreshold())
1791 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1792 if (dataQA->GetNTimeBins())
1793 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1794 if (dataQA->GetNPads())
1795 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1796 if (dataQA->GetTimePosition())
1797 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1799 prep.DumpToFile(file.Data());
1803 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1805 // Get the correction of the drift velocity using the laser tracks calbration
1808 // timeStamp - tim stamp in seconds
1809 // deltaT - integration period to calculate time0 offset
1810 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1811 // Note in case no data form both A and C side - the value from active side used
1812 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1813 TGraphErrors *grlaserA=0;
1814 TGraphErrors *grlaserC=0;
1815 Double_t vlaserA=0, vlaserC=0;
1816 if (!array) return 0;
1817 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1818 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1821 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1822 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1823 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1826 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1827 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1828 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1830 if (side==0) return vlaserA;
1831 if (side==1) return vlaserC;
1832 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1833 if (!grlaserA) return vlaserC;
1834 if (!grlaserC) return vlaserA;
1839 Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1841 // Get the correction of the drift velocity using the CE laser data
1842 // combining information from the CE, laser track calibration
1843 // and P/T calibration
1846 // timeStamp - tim stamp in seconds
1847 // deltaT - integration period to calculate time0 offset
1848 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1849 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
1850 if (!arrT) return 0;
1851 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1852 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1853 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1856 Double_t corrPTA = 0, corrPTC=0;
1857 Double_t ltime0A = 0, ltime0C=0;
1859 Double_t corrA=0, corrC=0;
1860 Double_t timeA=0, timeC=0;
1861 TGraph *graphA = (TGraph*)arrT->At(72);
1862 TGraph *graphC = (TGraph*)arrT->At(73);
1863 if (!graphA && !graphC) return 0.;
1864 if (graphA &&graphA->GetN()>0) {
1865 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1866 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1867 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1868 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1869 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1870 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1873 if (graphC&&graphC->GetN()>0){
1874 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1875 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1876 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1877 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1878 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1879 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1883 if (side ==0 ) return corrA;
1884 if (side ==1 ) return corrC;
1885 Double_t corrM= (corrA+corrC)*0.5;
1886 if (!graphA) corrM=corrC;
1887 if (!graphC) corrM=corrA;
1894 Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1896 // VERY obscure method - we need something in framework
1897 // Find the TPC runs with temperature OCDB entry
1898 // cache the start and end of the run
1900 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1901 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1902 if (!storage) return 0;
1903 TString path=storage->GetURI();
1907 if (path.Contains("local")){ // find the list if local system
1908 path.ReplaceAll("local://","");
1909 path+="TPC/Calib/Temperature";
1910 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
1912 runsT=gSystem->GetFromPipe(command);
1914 TObjArray *arr= runsT.Tokenize("r");
1917 TArrayI indexes(arr->GetEntries());
1918 TArrayI runs(arr->GetEntries());
1920 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1921 Int_t irunN = atoi(arr->At(irun)->GetName());
1922 if (irunN<startRun) continue;
1923 if (irunN>stopRun) continue;
1924 runs[naccept]=irunN;
1928 fRunsStart.Set(fRuns.fN);
1929 fRunsStop.Set(fRuns.fN);
1930 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1931 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
1934 AliCDBEntry * entry = 0;
1935 {for (Int_t irun=0;irun<fRuns.fN; irun++){
1936 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1937 if (!entry) continue;
1938 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1939 if (!tmpRun) continue;
1940 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1941 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1942 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1948 Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
1950 // binary search - find the run for given time stamp
1952 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
1953 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
1955 for (Int_t index=index0; index<=index1; index++){
1956 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
1958 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
1961 if (cindex<0) cindex =(index0+index1)/2;
1965 return fRuns[cindex];
1972 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
1974 // filter outlyer measurement
1975 // Only points around median +- sigmaCut filtered
1976 if (!graph) return 0;
1978 Int_t npoints0 = graph->GetN();
1981 Double_t *outx=new Double_t[npoints0];
1982 Double_t *outy=new Double_t[npoints0];
1985 if (npoints0<kMinPoints) return 0;
1986 for (Int_t iter=0; iter<3; iter++){
1988 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1989 if (graph->GetY()[ipoint]==0) continue;
1990 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
1991 outx[npoints] = graph->GetX()[ipoint];
1992 outy[npoints] = graph->GetY()[ipoint];
1995 if (npoints<=1) break;
1996 medianY =TMath::Median(npoints,outy);
1997 rmsY =TMath::RMS(npoints,outy);
2000 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2005 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2007 // filter outlyer measurement
2008 // Only points around median +- cut filtered
2009 if (!graph) return 0;
2011 Int_t npoints0 = graph->GetN();
2014 Double_t *outx=new Double_t[npoints0];
2015 Double_t *outy=new Double_t[npoints0];
2018 if (npoints0<kMinPoints) return 0;
2019 for (Int_t iter=0; iter<3; iter++){
2021 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2022 if (graph->GetY()[ipoint]==0) continue;
2023 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2024 outx[npoints] = graph->GetX()[ipoint];
2025 outy[npoints] = graph->GetY()[ipoint];
2028 if (npoints<=1) break;
2029 medianY =TMath::Median(npoints,outy);
2030 rmsY =TMath::RMS(npoints,outy);
2033 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2039 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
2041 // filter outlyer measurement
2042 // Only points with normalized errors median +- sigmaCut filtered
2044 Int_t kMinPoints=10;
2045 Int_t npoints0 = graph->GetN();
2047 Float_t medianErr=0, rmsErr=0;
2048 Double_t *outx=new Double_t[npoints0];
2049 Double_t *outy=new Double_t[npoints0];
2050 Double_t *erry=new Double_t[npoints0];
2051 Double_t *nerry=new Double_t[npoints0];
2052 Double_t *errx=new Double_t[npoints0];
2055 if (npoints0<kMinPoints) return 0;
2056 for (Int_t iter=0; iter<3; iter++){
2058 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2059 nerry[npoints] = graph->GetErrorY(ipoint);
2060 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2061 erry[npoints] = graph->GetErrorY(ipoint);
2062 outx[npoints] = graph->GetX()[ipoint];
2063 outy[npoints] = graph->GetY()[ipoint];
2064 errx[npoints] = graph->GetErrorY(ipoint);
2067 if (npoints==0) break;
2068 medianErr=TMath::Median(npoints,erry);
2069 medianY =TMath::Median(npoints,outy);
2070 rmsErr =TMath::RMS(npoints,erry);
2072 TGraphErrors *graphOut=0;
2073 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2082 void AliTPCcalibDButil::Sort(TGraph *graph){
2084 // sort array - neccessay for approx
2086 Int_t npoints = graph->GetN();
2087 Int_t *indexes=new Int_t[npoints];
2088 Double_t *outx=new Double_t[npoints];
2089 Double_t *outy=new Double_t[npoints];
2090 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2091 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2092 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2093 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2094 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2097 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2099 // smmoth graph - mean on the interval
2102 Int_t npoints = graph->GetN();
2103 Double_t *outy=new Double_t[npoints];
2105 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2106 Double_t lx=graph->GetX()[ipoint];
2107 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2108 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2109 if (index0<0) index0=0;
2110 if (index1>=npoints-1) index1=npoints-1;
2111 if ((index1-index0)>1){
2112 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2114 outy[ipoint]=graph->GetY()[ipoint];
2117 // TLinearFitter fitter(3,"pol2");
2118 // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2119 // Double_t lx=graph->GetX()[ipoint];
2120 // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2121 // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2122 // if (index0<0) index0=0;
2123 // if (index1>=npoints-1) index1=npoints-1;
2124 // fitter.ClearPoints();
2125 // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2126 // if ((index1-index0)>1){
2127 // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2129 // outy[ipoint]=graph->GetY()[ipoint];
2135 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2136 graph->GetY()[ipoint] = outy[ipoint];
2141 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
2143 // Use constant interpolation outside of range
2146 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2149 if (graph->GetN()<1){
2150 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2153 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2154 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2155 return graph->Eval( xref);
2158 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2160 // Filter DCS sensor information
2161 // ymin - minimal value
2163 // maxdy - maximal deirivative
2164 // sigmaCut - cut on values and derivative in terms of RMS distribution
2165 // Return value - accepted fraction
2169 // 0. Calculate median and rms of values in specified range
2170 // 1. Filter out outliers - median+-sigmaCut*rms
2171 // values replaced by median
2173 AliSplineFit * fit = sensor->GetFit();
2174 if (!fit) return 0.;
2175 Int_t nknots = fit->GetKnots();
2182 Double_t *yin0 = new Double_t[nknots];
2183 Double_t *yin1 = new Double_t[nknots];
2186 for (Int_t iknot=0; iknot< nknots; iknot++){
2187 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2188 yin0[naccept] = fit->GetY0()[iknot];
2189 yin1[naccept] = fit->GetY1()[iknot];
2190 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2200 Double_t medianY0=0, medianY1=0;
2201 Double_t rmsY0 =0, rmsY1=0;
2202 medianY0 = TMath::Median(naccept, yin0);
2203 medianY1 = TMath::Median(naccept, yin1);
2204 rmsY0 = TMath::RMS(naccept, yin0);
2205 rmsY1 = TMath::RMS(naccept, yin1);
2208 // 1. Filter out outliers - median+-sigmaCut*rms
2209 // values replaced by median
2210 // if replaced the derivative set to 0
2212 for (Int_t iknot=0; iknot< nknots; iknot++){
2214 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2215 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2216 if (nknots<2) fit->GetY1()[iknot]=0;
2217 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2219 fit->GetY0()[iknot]=medianY0;
2220 fit->GetY1()[iknot]=0;
2227 return Float_t(naccept)/Float_t(nknots);
2230 Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2232 // Filter temperature array
2233 // tempArray - array of temperatures -
2234 // ymin - minimal accepted temperature - default 15
2235 // ymax - maximal accepted temperature - default 22
2236 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2237 // return value - fraction of filtered sensors
2238 const Double_t kMaxDy=0.1;
2239 Int_t nsensors=tempArray->NumSensors();
2240 if (nsensors==0) return 0.;
2242 for (Int_t isensor=0; isensor<nsensors; isensor++){
2243 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2244 if (!sensor) continue;
2245 //printf("%d\n",isensor);
2246 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2247 if (sensor->GetFit()==0){
2249 tempArray->RemoveSensorNum(isensor);
2254 return Float_t(naccept)/Float_t(nsensors);
2258 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
2261 // Input parameters:
2262 // deltaT - smoothing window (in seconds)
2263 // cutAbs - max distance of the time info to the median (in time bins)
2264 // cutSigma - max distance (in the RMS)
2265 // pcstream - optional debug streamer to store original and filtered info
2266 // Hardwired parameters:
2267 // kMinPoints =10; // minimal number of points to define the CE
2268 // kMinSectors=12; // minimal number of sectors to define sideCE
2270 // 0. Filter almost emty graphs (kMinPoints=10)
2271 // 1. calculate median and RMS per side
2272 // 2. Filter graphs - in respect with side medians
2273 // - cutAbs and cutDelta used
2274 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2275 // 4. Calculate mean for A side and C side
2277 const Int_t kMinPoints =10; // minimal number of points to define the CE
2278 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2279 const Int_t kMinTime =400; // minimal arrival time of CE
2280 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2282 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2283 if (!cearray) return;
2288 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2289 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernPressure");
2290 if ( tempMapCE && cavernPressureCE){
2292 Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2293 FilterSensor(cavernPressureCE,960,1050,10, 5.);
2294 if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2296 // recalculate P/T correction map for time of the CE
2297 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2298 driftCalib->SetName("driftPTCE");
2299 driftCalib->SetTitle("driftPTCE");
2300 cearray->AddLast(driftCalib);
2304 // 0. Filter almost emty graphs
2307 for (Int_t i=0; i<72;i++){
2308 TGraph *graph= (TGraph*)arrT->At(i);
2309 if (!graph) continue;
2310 if (graph->GetN()<kMinPoints){
2312 delete graph; // delete empty graph
2315 if (tmin<0) tmin = graph->GetX()[0];
2316 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2318 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2319 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2322 // 1. calculate median and RMS per side
2324 TArrayF arrA(100000), arrC(100000);
2326 Double_t medianA=0, medianC=0;
2327 Double_t rmsA=0, rmsC=0;
2328 for (Int_t isec=0; isec<72;isec++){
2329 TGraph *graph= (TGraph*)arrT->At(isec);
2330 if (!graph) continue;
2331 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2332 if (graph->GetY()[ipoint]<kMinTime) continue;
2333 if (nA>=arrA.fN) arrA.Set(nA*2);
2334 if (nC>=arrC.fN) arrC.Set(nC*2);
2335 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2336 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2340 medianA=TMath::Median(nA,arrA.fArray);
2341 rmsA =TMath::RMS(nA,arrA.fArray);
2344 medianC=TMath::Median(nC,arrC.fArray);
2345 rmsC =TMath::RMS(nC,arrC.fArray);
2348 // 2. Filter graphs - in respect with side medians
2350 TArrayD vecX(100000), vecY(100000);
2351 for (Int_t isec=0; isec<72;isec++){
2352 TGraph *graph= (TGraph*)arrT->At(isec);
2353 if (!graph) continue;
2354 Double_t median = (isec%36<18) ? medianA: medianC;
2355 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2357 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2358 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2359 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2360 vecX[naccept]= graph->GetX()[ipoint];
2361 vecY[naccept]= graph->GetY()[ipoint];
2364 if (naccept<kMinPoints){
2365 arrT->AddAt(0,isec);
2366 delete graph; // delete empty graph
2369 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2371 arrT->AddAt(graph2,isec);
2374 // 3. Cut in respect wit the graph median
2376 for (Int_t i=0; i<72;i++){
2377 TGraph *graph= (TGraph*)arrT->At(i);
2378 if (!graph) continue;
2382 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2383 if (!graphTS0) continue;
2384 if (graphTS0->GetN()<kMinPoints) {
2390 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2392 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2394 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2395 (*pcstream)<<"filterCE"<<
2400 "graphTS0.="<<graphTS0<<
2401 "graphTS.="<<graphTS<<
2405 if (!graphTS) continue;
2406 arrT->AddAt(graphTS,i);
2410 // Recalculate the mean time A side C side
2412 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2413 Int_t meanPoints=(nA+nC)/72; // mean number of points
2414 for (Int_t itime=0; itime<200; itime++){
2416 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2417 for (Int_t i=0; i<72;i++){
2418 TGraph *graph= (TGraph*)arrT->At(i);
2419 if (!graph) continue;
2420 if (graph->GetN()<(meanPoints/4)) continue;
2421 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2422 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2426 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2427 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2428 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2429 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2432 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2433 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2435 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2436 (*pcstream)<<"filterAC"<<
2445 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2446 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2447 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2448 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2449 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2450 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2453 if (nA<kMinSectors) arrT->AddAt(0,72);
2454 else arrT->AddAt(graphTSA,72);
2455 if (nC<kMinSectors) arrT->AddAt(0,73);
2456 else arrT->AddAt(graphTSC,73);
2460 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
2462 // Filter Drift velocity measurement using the tracks
2463 // 0. remove outlyers - error based
2467 const Int_t kMinPoints=1; // minimal number of points to define value
2468 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2471 for (Int_t i=0; i<arrT->GetEntries();i++){
2472 TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2473 if (!graph) continue;
2474 if (graph->GetN()<kMinPoints){
2479 TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2481 delete graph; arrT->AddAt(0,i); continue;
2483 if (graph2->GetN()<1) {
2484 delete graph; arrT->AddAt(0,i); continue;
2486 graph2->SetName(graph->GetName());
2487 graph2->SetTitle(graph->GetTitle());
2488 arrT->AddAt(graph2,i);
2490 (*pcstream)<<"filterTracks"<<
2495 "graph2.="<<graph2<<
2506 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2509 // get laser time offset
2510 // median around timeStamp+-deltaT
2511 // QA - chi2 needed for later usage - to be added
2512 // - currently cut on error
2515 Double_t kMinDelay=0.01;
2516 Double_t kMinDelayErr=0.0001;
2518 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2519 if (!array) return 0;
2520 TGraphErrors *tlaser=0;
2522 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2523 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2525 if (!tlaser) return 0;
2526 Int_t npoints0= tlaser->GetN();
2527 if (npoints0==0) return 0;
2528 Double_t *xlaser = new Double_t[npoints0];
2529 Double_t *ylaser = new Double_t[npoints0];
2531 for (Int_t i=0;i<npoints0;i++){
2533 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2534 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2535 xlaser[npoints]=tlaser->GetX()[npoints];
2536 ylaser[npoints]=tlaser->GetY()[npoints];
2541 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2542 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2543 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2544 if (index0<0) index0=0;
2545 if (index1>=npoints-1) index1=npoints-1;
2546 if (index1-index0<kMinPoints) return 0;
2548 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2549 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2558 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector *pcstream){
2560 // Filter Goofie data
2561 // goofieArray - points will be filtered
2562 // deltaT - smmothing time window
2563 // cutSigma - outler sigma cut in rms
2564 // minVn, maxVd- range absolute cut for variable vd/pt
2567 // Ignore goofie if not enough points
2569 const Int_t kMinPoints = 3;
2572 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2573 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2574 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2575 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2576 if (!graphvd) return;
2577 if (graphvd->GetN()<kMinPoints){
2579 goofieArray->GetSensorNum(2)->SetGraph(0);
2583 // 1. Caluclate medians of critical variables
2589 Double_t medianpt=0;
2590 Double_t medianvd=0, sigmavd=0;
2591 Double_t medianan=0;
2592 Double_t medianaf=0;
2593 Int_t entries=graphvd->GetN();
2594 Double_t yvdn[10000];
2597 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2598 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2599 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2600 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2601 yvdn[nvd++]=graphvd->GetY()[ipoint];
2603 if (nvd<kMinPoints){
2605 goofieArray->GetSensorNum(2)->SetGraph(0);
2609 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2610 if (nuni>=kMinPoints){
2611 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2613 medianvd = TMath::Median(nvd, yvdn);
2616 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2617 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2618 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2619 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2620 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2621 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2629 // 2. Make outlyer graph
2632 TGraph graphOut(*graphvd);
2633 for (Int_t i=0; i<entries;i++){
2635 Bool_t isOut=kFALSE;
2636 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2637 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2639 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2641 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2642 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2643 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2644 graphOut.GetY()[i]= (isOut)?1:0;
2647 if (nOK<kMinPoints) {
2649 goofieArray->GetSensorNum(2)->SetGraph(0);
2653 // 3. Filter out outlyers - and smooth
2655 TVectorF vmedianArray(goofieArray->NumSensors());
2656 TVectorF vrmsArray(goofieArray->NumSensors());
2657 Double_t xnew[10000];
2658 Double_t ynew[10000];
2660 junk.SetOwner(kTRUE);
2664 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2666 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2667 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2669 if (!sensor) continue;
2670 graphOld = sensor->GetGraph();
2672 sensor->SetGraph(0);
2674 for (Int_t i=0;i<entries;i++){
2675 if (graphOut.GetY()[i]>0.5) continue;
2676 xnew[nused]=graphOld->GetX()[i];
2677 ynew[nused]=graphOld->GetY()[i];
2680 graphNew = new TGraph(nused,xnew,ynew);
2681 junk.AddLast(graphNew);
2682 junk.AddLast(graphOld);
2684 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2686 junk.AddLast(graphNew0);
2687 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2689 junk.AddLast(graphNew1);
2690 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2692 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2693 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2694 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2695 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2696 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2697 vmedianArray[isensor]=median;
2703 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2704 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2705 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2706 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2707 (*pcstream)<<"goofieA"<<
2708 Form("isOK_%d.=",isensor)<<isOK<<
2709 Form("s_%d.=",isensor)<<sensor<<
2710 Form("gr_%d.=",isensor)<<graphOld<<
2711 Form("gr0_%d.=",isensor)<<graphNew0<<
2712 Form("gr1_%d.=",isensor)<<graphNew1<<
2713 Form("gr2_%d.=",isensor)<<graphNew2;
2714 if (isOK) sensor->SetGraph(graphNew2);
2716 (*pcstream)<<"goofieA"<<
2717 "vmed.="<<&vmedianArray<<
2718 "vrms.="<<&vrmsArray<<
2720 junk.Delete(); // delete temoprary graphs