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>
33 #include <AliCDBStorage.h>
34 #include <AliDCSSensorArray.h>
35 #include <AliTPCSensorTempArray.h>
36 #include <AliDCSSensor.h>
37 #include "AliTPCcalibDB.h"
38 #include "AliTPCCalPad.h"
39 #include "AliTPCCalROC.h"
40 #include "AliTPCROC.h"
41 #include "AliTPCmapper.h"
42 #include "AliTPCParam.h"
43 #include "AliTPCCalibRaw.h"
44 #include "TGraphErrors.h"
46 #include "AliTPCcalibDButil.h"
47 #include "AliTPCPreprocessorOnline.h"
48 #include "AliTPCCalibVdrift.h"
49 #include "AliMathBase.h"
51 ClassImp(AliTPCcalibDButil)
52 AliTPCcalibDButil::AliTPCcalibDButil() :
60 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
71 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
78 fMapper(new AliTPCmapper(0x0)),
82 fPulTmaxLimitAbs(1.5),
85 fRuns(0), // run list with OCDB info
86 fRunsStart(0), // start time for given run
87 fRunsStop(0) // stop time for given run
93 //_____________________________________________________________________________________
94 AliTPCcalibDButil::~AliTPCcalibDButil()
99 delete fPulserOutlier;
100 delete fRefPulserOutlier;
102 if (fRefPadNoise) delete fRefPadNoise;
103 if (fRefPedestals) delete fRefPedestals;
104 if (fRefPulserTmean) delete fRefPulserTmean;
105 if (fRefPulserTrms) delete fRefPulserTrms;
106 if (fRefPulserQmean) delete fRefPulserQmean;
107 if (fRefCETmean) delete fRefCETmean;
108 if (fRefCETrms) delete fRefCETrms;
109 if (fRefCEQmean) delete fRefCEQmean;
110 if (fRefALTROMasked) delete fRefALTROMasked;
111 if (fRefCalibRaw) delete fRefCalibRaw;
114 //_____________________________________________________________________________________
115 void AliTPCcalibDButil::UpdateFromCalibDB()
118 // Update pointers from calibDB
120 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
121 fPadNoise=fCalibDB->GetPadNoise();
122 fPedestals=fCalibDB->GetPedestals();
123 fPulserTmean=fCalibDB->GetPulserTmean();
124 fPulserTrms=fCalibDB->GetPulserTrms();
125 fPulserQmean=fCalibDB->GetPulserQmean();
126 fCETmean=fCalibDB->GetCETmean();
127 fCETrms=fCalibDB->GetCETrms();
128 fCEQmean=fCalibDB->GetCEQmean();
129 fALTROMasked=fCalibDB->GetALTROMasked();
130 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
131 fCalibRaw=fCalibDB->GetCalibRaw();
132 UpdatePulserOutlierMap();
134 //_____________________________________________________________________________________
135 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
136 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
139 // Process the CE data for this run
140 // the return TVectorD arrays contian the results of the fit
141 // noutliersCE contains the number of pads marked as outliers,
142 // not including masked and edge pads
145 //retrieve CE and ALTRO data
147 TString fitString(fitFormula);
148 fitString.ReplaceAll("++","#");
149 Int_t ndim=fitString.CountChar('#')+2;
150 fitResultsA.ResizeTo(ndim);
151 fitResultsC.ResizeTo(ndim);
160 if (outCE) out=outCE;
161 else out=new AliTPCCalPad("outCE","outCE");
162 AliTPCCalROC *rocMasked=0x0;
163 //loop over all channels
164 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
165 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
166 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
167 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
169 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
173 //add time offset to IROCs
174 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
175 rocData->Add(fIrocTimeOffset);
177 UInt_t nrows=rocData->GetNrows();
178 for (UInt_t irow=0;irow<nrows;++irow){
179 UInt_t npads=rocData->GetNPads(irow);
180 for (UInt_t ipad=0;ipad<npads;++ipad){
181 rocOut->SetValue(irow,ipad,0);
182 //exclude masked pads
183 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
184 rocOut->SetValue(irow,ipad,1);
187 //exclude first two rows in IROC and last two rows in OROC
189 if (irow<2) rocOut->SetValue(irow,ipad,1);
191 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
194 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
195 Float_t valTmean=rocData->GetValue(irow,ipad);
196 //exclude values that are exactly 0
198 rocOut->SetValue(irow,ipad,1);
201 // exclude channels with too large variations
202 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
203 rocOut->SetValue(irow,ipad,1);
211 Float_t chi2Af,chi2Cf;
212 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
215 if (!outCE) delete out;
217 //_____________________________________________________________________________________
218 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
219 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
220 Float_t &driftTimeA, Float_t &driftTimeC )
223 // Calculate statistical information from the CE graphs for drift time and charge
227 vecTEntries.ResizeTo(72);
228 vecTMean.ResizeTo(72);
229 vecTRMS.ResizeTo(72);
230 vecTMedian.ResizeTo(72);
231 vecQEntries.ResizeTo(72);
232 vecQMean.ResizeTo(72);
233 vecQRMS.ResizeTo(72);
234 vecQMedian.ResizeTo(72);
245 TObjArray *arrT=fCalibDB->GetCErocTtime();
246 TObjArray *arrQ=fCalibDB->GetCErocQtime();
248 for (Int_t isec=0;isec<74;++isec){
249 TGraph *gr=(TGraph*)arrT->At(isec);
252 Int_t npoints = gr->GetN();
253 values.ResizeTo(npoints);
255 //skip first points, theres always some problems with finding the CE position
256 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
257 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
258 values[nused]=gr->GetY()[ipoint];
263 if (isec<72) vecTEntries[isec]= nused;
266 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
267 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
268 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
269 } else if (isec==72){
270 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
271 } else if (isec==73){
272 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
278 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
279 TGraph *gr=(TGraph*)arrQ->At(isec);
282 Int_t npoints = gr->GetN();
283 values.ResizeTo(npoints);
285 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
286 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
287 values[nused]=gr->GetY()[ipoint];
292 vecQEntries[isec]= nused;
294 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
295 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
296 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
302 //_____________________________________________________________________________________
303 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
304 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
305 Int_t &nonMaskedZero)
308 // process noise data
309 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
310 // OROCs small pads [2] and OROCs large pads [3]
311 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
312 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
315 //set proper size and reset
316 const UInt_t infoSize=4;
317 vNoiseMean.ResizeTo(infoSize);
318 vNoiseMeanSenRegions.ResizeTo(infoSize);
319 vNoiseRMS.ResizeTo(infoSize);
320 vNoiseRMSSenRegions.ResizeTo(infoSize);
322 vNoiseMeanSenRegions.Zero();
324 vNoiseRMSSenRegions.Zero();
327 TVectorD c(infoSize);
328 TVectorD cs(infoSize);
332 //retrieve noise and ALTRO data
333 if (!fPadNoise) return;
334 AliTPCCalROC *rocMasked=0x0;
335 //create IROC, OROC1, OROC2 and sensitive region masks
336 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
337 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
338 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
339 UInt_t nrows=noiseROC->GetNrows();
340 for (UInt_t irow=0;irow<nrows;++irow){
341 UInt_t npads=noiseROC->GetNPads(irow);
342 for (UInt_t ipad=0;ipad<npads;++ipad){
343 //don't use masked channels;
344 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
345 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
352 if ( !(noiseVal<10000000) ){
353 printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
356 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
357 Int_t masksen=1; // sensitive pards are not masked (0)
358 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
359 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
361 if (irow>19&&irow<46){
362 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
365 vNoiseMean[type]+=noiseVal;
366 vNoiseRMS[type]+=noiseVal*noiseVal;
369 vNoiseMeanSenRegions[type]+=noiseVal;
370 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
375 //define sensive regions
376 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
378 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
379 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
381 if ((Int_t)irow<par.GetNRowUp1()){
384 vNoiseMean[type]+=noiseVal;
385 vNoiseRMS[type]+=noiseVal*noiseVal;
388 vNoiseMeanSenRegions[type]+=noiseVal;
389 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
395 vNoiseMean[type]+=noiseVal;
396 vNoiseRMS[type]+=noiseVal*noiseVal;
399 vNoiseMeanSenRegions[type]+=noiseVal;
400 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
407 vNoiseMean[type]+=noiseVal;
408 vNoiseRMS[type]+=noiseVal*noiseVal;
411 vNoiseMeanSenRegions[type]+=noiseVal;
412 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
417 }//end loop sectors (rocs)
419 //calculate mean and RMS
420 const Double_t verySmall=0.0000000001;
421 for (UInt_t i=0;i<infoSize;++i){
428 // printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
429 mean=vNoiseMean[i]/c[i];
431 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
436 if (cs[i]>verySmall){
437 meanSen=vNoiseMeanSenRegions[i]/cs[i];
438 rmsSen=vNoiseRMSSenRegions[i];
439 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
441 vNoiseMeanSenRegions[i]=meanSen;
442 vNoiseRMSSenRegions[i]=rmsSen;
446 //_____________________________________________________________________________________
447 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
450 // Process the Pulser information
451 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
454 const UInt_t infoSize=4;
455 //reset counters to error number
456 vMeanTime.ResizeTo(infoSize);
459 TVectorD c(infoSize);
460 //retrieve pulser and ALTRO data
461 if (!fPulserTmean) return;
464 AliTPCCalROC *rocOut=0x0;
465 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
466 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
467 if (!tmeanROC) continue;
468 rocOut=fPulserOutlier->GetCalROC(isec);
469 UInt_t nchannels=tmeanROC->GetNchannels();
470 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
471 if (rocOut && rocOut->GetValue(ichannel)) continue;
472 Float_t val=tmeanROC->GetValue(ichannel);
474 vMeanTime[type]+=val;
479 for (UInt_t itype=0; itype<infoSize; ++itype){
480 if (c[itype]>0) vMeanTime[itype]/=c[itype];
481 else vMeanTime[itype]=0;
484 //_____________________________________________________________________________________
485 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
488 // Get Values from ALTRO configuration data
491 if (!fALTROMasked) return;
493 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
494 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
495 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
496 if (rocMasked->GetValue(ichannel)) ++nMasked;
500 //_____________________________________________________________________________________
501 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
504 // Proces Goofie values, return statistical information of the currently set goofieArray
505 // The meaning of the entries are given below
507 1 TPC_ANODE_I_A00_STAT
509 3 TPC_DVM_DriftVelocity
514 8 TPC_DVM_NumberOfSparks
515 9 TPC_DVM_PeakAreaFar
516 10 TPC_DVM_PeakAreaNear
517 11 TPC_DVM_PeakPosFar
518 12 TPC_DVM_PeakPosNear
524 18 TPC_DVM_TemperatureS1
528 vecEntries.ResizeTo(nsensors);
529 vecMedian.ResizeTo(nsensors);
530 vecMean.ResizeTo(nsensors);
531 vecRMS.ResizeTo(nsensors);
538 Double_t kEpsilon=0.0000000001;
539 Double_t kBig=100000000000.;
540 Int_t nsensors = fGoofieArray->NumSensors();
541 vecEntries.ResizeTo(nsensors);
542 vecMedian.ResizeTo(nsensors);
543 vecMean.ResizeTo(nsensors);
544 vecRMS.ResizeTo(nsensors);
546 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
547 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
548 if (gsensor && gsensor->GetGraph()){
549 Int_t npoints = gsensor->GetGraph()->GetN();
551 values.ResizeTo(npoints);
553 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
554 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
555 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
556 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
561 vecEntries[isensor]= nused;
563 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
564 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
565 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
570 //_____________________________________________________________________________________
571 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
574 // check the variations of the pedestal data to the reference pedestal data
575 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
578 TVectorF vThres(npar); //thresholds
579 Int_t nActive=0; //number of active channels
581 //reset and set thresholds
582 pedestalDeviations.ResizeTo(npar);
583 for (Int_t i=0;i<npar;++i){
584 pedestalDeviations.GetMatrixArray()[i]=0;
585 vThres.GetMatrixArray()[i]=(i+1)*.5;
587 //check all needed data is available
588 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
589 //loop over all channels
590 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
591 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
592 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
593 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
594 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
595 UInt_t nrows=mROC->GetNrows();
596 for (UInt_t irow=0;irow<nrows;++irow){
597 UInt_t npads=mROC->GetNPads(irow);
598 for (UInt_t ipad=0;ipad<npads;++ipad){
599 //don't use masked channels;
600 if (mROC ->GetValue(irow,ipad)) continue;
601 if (mRefROC->GetValue(irow,ipad)) continue;
602 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
603 for (Int_t i=0;i<npar;++i){
604 if (deviation>vThres[i])
605 ++pedestalDeviations.GetMatrixArray()[i];
612 for (Int_t i=0;i<npar;++i){
613 pedestalDeviations.GetMatrixArray()[i]/=nActive;
617 //_____________________________________________________________________________________
618 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
621 // check the variations of the noise data to the reference noise data
622 // thresholds are 5, 10, 15 and 20 percent respectively.
625 TVectorF vThres(npar); //thresholds
626 Int_t nActive=0; //number of active channels
628 //reset and set thresholds
629 noiseDeviations.ResizeTo(npar);
630 for (Int_t i=0;i<npar;++i){
631 noiseDeviations.GetMatrixArray()[i]=0;
632 vThres.GetMatrixArray()[i]=(i+1)*.05;
634 //check all needed data is available
635 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
636 //loop over all channels
637 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
638 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
639 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
640 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
641 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
642 UInt_t nrows=mROC->GetNrows();
643 for (UInt_t irow=0;irow<nrows;++irow){
644 UInt_t npads=mROC->GetNPads(irow);
645 for (UInt_t ipad=0;ipad<npads;++ipad){
646 //don't use masked channels;
647 if (mROC ->GetValue(irow,ipad)) continue;
648 if (mRefROC->GetValue(irow,ipad)) continue;
649 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
650 for (Int_t i=0;i<npar;++i){
651 if (deviation>vThres[i])
652 ++noiseDeviations.GetMatrixArray()[i];
659 for (Int_t i=0;i<npar;++i){
660 noiseDeviations.GetMatrixArray()[i]/=nActive;
664 //_____________________________________________________________________________________
665 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
666 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
669 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
670 // thresholds are .5, 1, 5 and 10 percent respectively.
674 TVectorF vThres(npar); //thresholds
675 Int_t nActive=0; //number of active channels
677 //reset and set thresholds
678 pulserQdeviations.ResizeTo(npar);
679 for (Int_t i=0;i<npar;++i){
680 pulserQdeviations.GetMatrixArray()[i]=0;
685 vThres.GetMatrixArray()[0]=.005;
686 vThres.GetMatrixArray()[1]=.01;
687 vThres.GetMatrixArray()[2]=.05;
688 vThres.GetMatrixArray()[3]=.1;
689 //check all needed data is available
690 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
692 UpdateRefPulserOutlierMap();
693 //loop over all channels
694 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
695 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
696 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
697 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
698 // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
699 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
700 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
701 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
702 Float_t pt_mean=ptROC->GetMean(oROC);
703 UInt_t nrows=mROC->GetNrows();
704 for (UInt_t irow=0;irow<nrows;++irow){
705 UInt_t npads=mROC->GetNPads(irow);
706 for (UInt_t ipad=0;ipad<npads;++ipad){
707 //don't use masked channels;
708 if (mROC ->GetValue(irow,ipad)) continue;
709 if (mRefROC->GetValue(irow,ipad)) continue;
710 //don't user edge pads
711 if (ipad==0||ipad==npads-1) continue;
713 Float_t pq=pqROC->GetValue(irow,ipad);
714 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
715 Float_t pt=ptROC->GetValue(irow,ipad);
716 // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
718 Float_t deviation=TMath::Abs(pq/pqRef-1);
719 for (Int_t i=0;i<npar;++i){
720 if (deviation>vThres[i])
721 ++pulserQdeviations.GetMatrixArray()[i];
723 if (pqRef>11&&pq<11) ++npadsOffAdd;
726 if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
732 for (Int_t i=0;i<npar;++i){
733 pulserQdeviations.GetMatrixArray()[i]/=nActive;
738 //_____________________________________________________________________________________
739 void AliTPCcalibDButil::UpdatePulserOutlierMap()
744 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
746 //_____________________________________________________________________________________
747 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
752 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
754 //_____________________________________________________________________________________
755 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
758 // Create a map that contains outliers from the Pulser calibration data.
759 // The outliers include masked channels, edge pads and pads with
760 // too large timing and charge variations.
761 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
762 // those do not contain masked and edge pads
766 pulOut->Multiply(0.);
770 AliTPCCalROC *rocMasked=0x0;
774 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
775 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
776 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
777 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
778 if (!tmeanROC||!qmeanROC) {
779 //reset outliers in this ROC
780 outROC->Multiply(0.);
783 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
785 // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
786 // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
787 UInt_t nrows=tmeanROC->GetNrows();
788 for (UInt_t irow=0;irow<nrows;++irow){
789 UInt_t npads=tmeanROC->GetNPads(irow);
790 for (UInt_t ipad=0;ipad<npads;++ipad){
791 Int_t outlier=0,masked=0;
792 Float_t q=qmeanROC->GetValue(irow,ipad);
793 Float_t t=tmeanROC->GetValue(irow,ipad);
794 //masked channels are outliers
795 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
796 //edge pads are outliers
797 if (ipad==0||ipad==npads-1) masked=1;
798 //channels with too large charge or timing deviation from the meadian are outliers
799 // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
800 if (q<fPulQminLimit && !masked) outlier=1;
802 if ( !(q<10000000) || !(t<10000000)) outlier=1;
803 outROC->SetValue(irow,ipad,outlier+masked);
804 fNpulserOutliers+=outlier;
809 //_____________________________________________________________________________________
810 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
813 // Create pad time0 object from pulser and/or CE data, depending on the selected model
814 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
815 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
816 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
818 // In case model 2 is invoked - gy arival time gradient is also returned
822 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
823 // decide between different models
824 if (model==0||model==1){
826 if (model==1) ProcessPulser(vMean);
827 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
828 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
829 if (!rocPulTmean) continue;
830 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
831 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
832 Float_t mean=rocPulTmean->GetMean(rocOut);
833 //treat case where a whole partition is masked
834 if (mean==0) mean=rocPulTmean->GetMean();
839 UInt_t nrows=rocTime0->GetNrows();
840 for (UInt_t irow=0;irow<nrows;++irow){
841 UInt_t npads=rocTime0->GetNPads(irow);
842 for (UInt_t ipad=0;ipad<npads;++ipad){
843 Float_t time=rocPulTmean->GetValue(irow,ipad);
844 //in case of an outlier pad use the mean of the altro values.
845 //This should be the most precise guess in that case.
846 if (rocOut->GetValue(irow,ipad)) {
847 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
848 if (time==0) time=mean;
850 Float_t val=time-mean;
851 rocTime0->SetValue(irow,ipad,val);
855 } else if (model==2){
856 Double_t pgya,pgyc,pchi2a,pchi2c;
857 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
858 fCETmean->Add(padPulser,-1.);
860 AliTPCCalPad outCE("outCE","outCE");
862 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
863 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
864 // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
865 if (!padFit) { delete padPulser; return 0;}
868 fCETmean->Add(padPulser,1.);
869 padTime0->Add(fCETmean);
870 padTime0->Add(padFit,-1);
875 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
876 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
877 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
878 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
879 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
880 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
881 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
882 Float_t mean=rocPulTmean->GetMean(rocOutPul);
883 if (mean==0) mean=rocPulTmean->GetMean();
884 UInt_t nrows=rocTime0->GetNrows();
885 for (UInt_t irow=0;irow<nrows;++irow){
886 UInt_t npads=rocTime0->GetNPads(irow);
887 for (UInt_t ipad=0;ipad<npads;++ipad){
888 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
889 if (rocOutCE->GetValue(irow,ipad)){
890 Float_t valOut=rocCEfit->GetValue(irow,ipad);
891 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
892 rocTime0->SetValue(irow,ipad,valOut);
900 Double_t median = padTime0->GetMedian();
901 padTime0->Add(-median); // normalize to median
904 //_____________________________________________________________________________________
905 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
907 if (roc==0) return 0.;
908 const Int_t sector=roc->GetSector();
909 AliTPCROC *tpcRoc=AliTPCROC::Instance();
910 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
914 //loop over a small range around the requested pad (+-10 rows/pads)
915 for (Int_t irow=row-10;irow<row+10;++irow){
916 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
917 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
918 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
919 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
920 if (altroRoc!=altroCurr) continue;
921 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
922 Float_t val=roc->GetValue(irow,ipad);
930 //_____________________________________________________________________________________
931 void AliTPCcalibDButil::SetRefFile(const char* filename)
934 // load cal pad objects form the reference file
936 TDirectory *currDir=gDirectory;
938 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
939 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
941 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
942 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
943 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
945 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
946 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
947 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
949 // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
950 // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
951 // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
952 // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
953 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
961 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
963 // Author: marian.ivanov@cern.ch
965 // Create outlier map for CE study
967 // Return value - outlyer map
968 // noutlyersCE - number of outlyers
969 // minSignal - minimal total Q signal
970 // cutRMSMin - minimal width of the signal in respect to the median
971 // cutRMSMax - maximal width of the signal in respect to the median
972 // cutMaxDistT - maximal deviation from time median per chamber
974 // Outlyers criteria:
975 // 0. Exclude masked pads
976 // 1. Exclude first two rows in IROC and last two rows in OROC
977 // 2. Exclude edge pads
978 // 3. Exclude channels with too large variations
979 // 4. Exclude pads with too small signal
980 // 5. Exclude signal with outlyers RMS
981 // 6. Exclude channels to far from the chamber median
984 AliTPCCalPad *out=ceOut;
985 if (!out) out= new AliTPCCalPad("outCE","outCE");
986 AliTPCCalROC *rocMasked=0x0;
987 if (!fCETmean) return 0;
988 if (!fCETrms) return 0;
989 if (!fCEQmean) return 0;
991 //loop over all channels
993 Double_t rmsMedian = fCETrms->GetMedian();
994 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
995 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
996 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
997 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
998 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
999 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1000 Double_t trocMedian = rocData->GetMedian();
1003 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1009 UInt_t nrows=rocData->GetNrows();
1010 for (UInt_t irow=0;irow<nrows;++irow){
1011 UInt_t npads=rocData->GetNPads(irow);
1012 for (UInt_t ipad=0;ipad<npads;++ipad){
1013 rocOut->SetValue(irow,ipad,0);
1014 Float_t valTmean=rocData->GetValue(irow,ipad);
1015 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1016 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1017 //0. exclude masked pads
1018 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1019 rocOut->SetValue(irow,ipad,1);
1022 //1. exclude first two rows in IROC and last two rows in OROC
1024 if (irow<2) rocOut->SetValue(irow,ipad,1);
1026 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1028 //2. exclude edge pads
1029 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1030 //exclude values that are exactly 0
1032 rocOut->SetValue(irow,ipad,1);
1035 //3. exclude channels with too large variations
1036 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1037 rocOut->SetValue(irow,ipad,1);
1041 //4. exclude channels with too small signal
1042 if (valQmean<minSignal) {
1043 rocOut->SetValue(irow,ipad,1);
1047 //5. exclude channels with too small rms
1048 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1049 rocOut->SetValue(irow,ipad,1);
1053 //6. exclude channels to far from the chamber median
1054 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1055 rocOut->SetValue(irow,ipad,1);
1066 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1068 // Author: marian.ivanov@cern.ch
1070 // Create outlier map for Pulser
1072 // Return value - outlyer map
1073 // noutlyersPulser - number of outlyers
1074 // cutTime - absolute cut - distance to the median of chamber
1075 // cutnRMSQ - nsigma cut from median q distribution per chamber
1076 // cutnRMSrms - nsigma cut from median rms distribution
1077 // Outlyers criteria:
1078 // 0. Exclude masked pads
1079 // 1. Exclude time outlyers (default 3 time bins)
1080 // 2. Exclude q outlyers (default 5 sigma)
1081 // 3. Exclude rms outlyers (default 5 sigma)
1083 AliTPCCalPad *out=pulserOut;
1084 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1085 AliTPCCalROC *rocMasked=0x0;
1086 if (!fPulserTmean) return 0;
1087 if (!fPulserTrms) return 0;
1088 if (!fPulserQmean) return 0;
1090 //loop over all channels
1092 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1093 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1094 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1095 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1096 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1097 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1099 Double_t rocMedianT = rocData->GetMedian();
1100 Double_t rocMedianQ = rocPulserQ->GetMedian();
1101 Double_t rocRMSQ = rocPulserQ->GetRMS();
1102 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1103 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1104 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1105 rocOut->SetValue(ichannel,0);
1106 Float_t valTmean=rocData->GetValue(ichannel);
1107 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1108 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1110 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1111 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1112 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1113 rocOut->SetValue(ichannel,isOut);
1114 if (isOut) noutliersPulser++;
1121 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1123 // Author : Marian Ivanov
1124 // Create pad time0 correction map using information from the CE and from pulser
1127 // Return PadTime0 to be used for time0 relative alignment
1128 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1129 // using $ALICE_ROOT/TPC/script/gui application
1131 // fitResultsA - fitParameters A side
1132 // fitResultsC - fitParameters C side
1133 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1134 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1138 // 1. Find outlier map for CE
1139 // 2. Find outlier map for Pulser
1140 // 3. Replace outlier by median at given sector (median without outliers)
1141 // 4. Substract from the CE data pulser
1142 // 5. Fit the CE with formula
1143 // 5.1) (IROC-OROC) offset
1147 // 5.5) (IROC-OROC)*(lx-xmid)
1149 // 6. Substract gy fit dependence from the CE data
1150 // 7. Add pulser back to CE data
1151 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1152 // 9. return CE data
1154 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1155 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1158 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)";
1159 // output for fit formula
1160 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)";
1161 // gy part of formula
1162 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)";
1165 if (!fCETmean) return 0;
1166 Double_t pgya,pgyc,pchi2a,pchi2c;
1167 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1168 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1170 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1171 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1172 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1173 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1174 padPulser->SetName("padPulser");
1175 padPulserOut->SetName("padPulserOut");
1176 padCE->SetName("padCE");
1177 padCEIn->SetName("padCEIn");
1178 padCEOut->SetName("padCEOut");
1179 padOut->SetName("padOut");
1182 // make combined outlyers map
1183 // and replace outlyers in maps with median for chamber
1185 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1186 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1187 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1188 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1189 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1190 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1191 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1192 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1193 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1194 if (rocPulserOut->GetValue(ichannel)>0) {
1195 rocPulser->SetValue(ichannel,pulserMedian);
1196 rocOut->SetValue(ichannel,1);
1198 if (rocCEOut->GetValue(ichannel)>0) {
1199 rocCE->SetValue(ichannel,ceMedian);
1200 rocOut->SetValue(ichannel,1);
1205 // remove pulser time 0
1207 padCE->Add(padPulser,-1);
1212 Float_t chi2Af,chi2Cf;
1213 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1217 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1218 padCEFitGY->SetName("padCEFitGy");
1220 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1221 padCEFit->SetName("padCEFit");
1223 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1224 padCEDiff->SetName("padCEDiff");
1225 padCEDiff->Add(padCEFit,-1.);
1228 padCE->Add(padCEFitGY,-1.);
1230 padCE->Add(padPulser,1.);
1231 Double_t padmedian = padCE->GetMedian();
1232 padCE->Add(-padmedian); // normalize to median
1234 // Replace outliers by fit value - median of diff per given chamber -GY fit
1236 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1237 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1238 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1239 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1240 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1241 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1243 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1244 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1245 if (rocOut->GetValue(ichannel)==0) continue;
1246 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1247 rocCE->SetValue(ichannel,value);
1253 //dump to the file - result can be visualized
1254 AliTPCPreprocessorOnline preprocesor;
1255 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1256 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1257 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1258 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1260 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1261 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1263 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1264 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1265 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1266 preprocesor.DumpToFile(dumpfile);
1269 delete padPulserOut;
1282 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1284 // find the closest point to xref in x direction
1285 // return dx and value
1287 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1288 if (index<0) index=0;
1289 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1290 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1291 dx = xref-graph->GetX()[index];
1292 y = graph->GetY()[index];
1297 Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1299 // Get the correction of the trigger offset
1300 // combining information from the laser track calibration
1301 // and from cosmic calibration
1304 // timeStamp - tim stamp in seconds
1305 // deltaT - integration period to calculate offset
1306 // deltaTLaser -max validity of laser data
1307 // valType - 0 - median, 1- mean
1309 // Integration vaues are just recomendation - if not possible to get points
1310 // automatically increase the validity by factor 2
1311 // (recursive algorithm until one month of data taking)
1314 const Float_t kLaserCut=0.0005;
1315 const Int_t kMaxPeriod=3600*24*30*3; // 3 month max
1316 const Int_t kMinPoints=20;
1318 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1320 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1322 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1323 if (!array) return 0;
1325 TGraphErrors *laserA[3]={0,0,0};
1326 TGraphErrors *laserC[3]={0,0,0};
1327 TGraphErrors *cosmicAll=0;
1328 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1329 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1330 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1333 if (!cosmicAll) return 0;
1334 Int_t nmeasC=cosmicAll->GetN();
1335 Float_t *tdelta = new Float_t[nmeasC];
1337 for (Int_t i=0;i<nmeasC;i++){
1338 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1339 Float_t ccosmic=cosmicAll->GetY()[i];
1340 Double_t yA=0,yC=0,dA=0,dC=0;
1341 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1342 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1343 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1344 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1346 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1348 if (TMath::Abs(yA-yC)<kLaserCut) {
1351 if (i%2==0) claser=yA;
1352 if (i%2==1) claser=yC;
1354 tdelta[nused]=ccosmic-claser;
1357 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1358 Double_t median = TMath::Median(nused,tdelta);
1359 Double_t mean = TMath::Mean(nused,tdelta);
1361 return (valType==0) ? median:mean;
1364 Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1366 // Get the correction of the drift velocity
1367 // combining information from the laser track calibration
1368 // and from cosmic calibration
1370 // dist - return value - distance to closest point in graph
1372 // timeStamp - tim stamp in seconds
1373 // deltaT - integration period to calculate time0 offset
1374 // deltaTLaser -max validity of laser data
1375 // valType - 0 - median, 1- mean
1377 // Integration vaues are just recomendation - if not possible to get points
1378 // automatically increase the validity by factor 2
1379 // (recursive algorithm until one month of data taking)
1383 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1385 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1387 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1388 if (!array) return 0;
1389 TGraphErrors *cosmicAll=0;
1390 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1391 if (!cosmicAll) return 0;
1393 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1395 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1396 Double_t vcosmic= AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1397 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1398 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1405 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1406 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1407 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1409 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1410 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1411 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1412 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1414 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1415 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1421 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1423 // Get the correction of the drift velocity using the laser tracks calbration
1426 // timeStamp - tim stamp in seconds
1427 // deltaT - integration period to calculate time0 offset
1428 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1429 // Note in case no data form both A and C side - the value from active side used
1430 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1431 TGraphErrors *grlaserA=0;
1432 TGraphErrors *grlaserC=0;
1433 Double_t vlaserA=0, vlaserC=0;
1434 if (!array) return 0;
1435 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1436 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1439 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1440 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1441 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1444 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1445 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1446 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1448 if (side==0) return vlaserA;
1449 if (side==1) return vlaserC;
1450 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1451 if (!grlaserA) return vlaserC;
1452 if (!grlaserC) return vlaserA;
1457 Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1459 // Get the correction of the drift velocity using the CE laser data
1460 // combining information from the CE, laser track calibration
1461 // and P/T calibration
1464 // timeStamp - tim stamp in seconds
1465 // deltaT - integration period to calculate time0 offset
1466 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1467 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
1468 if (!arrT) return 0;
1469 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1470 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1471 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1474 Double_t corrPTA = 0, corrPTC=0;
1475 Double_t ltime0A = 0, ltime0C=0;
1477 Double_t corrA=0, corrC=0;
1478 Double_t timeA=0, timeC=0;
1479 TGraph *graphA = (TGraph*)arrT->At(72);
1480 TGraph *graphC = (TGraph*)arrT->At(73);
1481 if (!graphA && !graphC) return 0.;
1482 if (graphA &&graphA->GetN()>0) {
1483 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1484 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1485 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1486 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1487 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1488 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1491 if (graphC&&graphC->GetN()>0){
1492 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1493 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1494 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1495 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1496 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1497 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1501 if (side ==0 ) return corrA;
1502 if (side ==1 ) return corrC;
1503 Double_t corrM= (corrA+corrC)*0.5;
1504 if (!graphA) corrM=corrC;
1505 if (!graphC) corrM=corrA;
1512 Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1514 // VERY obscure method - we need something in framework
1515 // Find the TPC runs with temperature OCDB entry
1516 // cache the start and end of the run
1518 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1519 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1520 if (!storage) return 0;
1521 TString path=storage->GetURI();
1525 if (path.Contains("local")){ // find the list if local system
1526 path.ReplaceAll("local://","");
1527 path+="TPC/Calib/Temperature";
1528 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
1530 runsT=gSystem->GetFromPipe(command);
1532 TObjArray *arr= runsT.Tokenize("r");
1535 TArrayI indexes(arr->GetEntries());
1536 TArrayI runs(arr->GetEntries());
1538 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1539 Int_t irunN = atoi(arr->At(irun)->GetName());
1540 if (irunN<startRun) continue;
1541 if (irunN>stopRun) continue;
1542 runs[naccept]=irunN;
1546 fRunsStart.Set(fRuns.fN);
1547 fRunsStop.Set(fRuns.fN);
1548 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1549 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
1552 AliCDBEntry * entry = 0;
1553 {for (Int_t irun=0;irun<fRuns.fN; irun++){
1554 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1555 if (!entry) continue;
1556 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1557 if (!tmpRun) continue;
1558 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1559 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1560 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1566 Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
1568 // binary search - find the run for given time stamp
1570 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
1571 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
1573 for (Int_t index=index0; index<=index1; index++){
1574 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
1576 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
1579 if (cindex<0) cindex =(index0+index1)/2;
1583 return fRuns[cindex];
1590 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
1592 // filter outlyer measurement
1593 // Only points around median +- sigmaCut filtered
1594 if (!graph) return 0;
1596 Int_t npoints0 = graph->GetN();
1599 Double_t *outx=new Double_t[npoints0];
1600 Double_t *outy=new Double_t[npoints0];
1603 if (npoints0<kMinPoints) return 0;
1604 for (Int_t iter=0; iter<3; iter++){
1606 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1607 if (graph->GetY()[ipoint]==0) continue;
1608 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
1609 outx[npoints] = graph->GetX()[ipoint];
1610 outy[npoints] = graph->GetY()[ipoint];
1613 if (npoints<=1) break;
1614 medianY =TMath::Median(npoints,outy);
1615 rmsY =TMath::RMS(npoints,outy);
1618 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1623 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
1625 // filter outlyer measurement
1626 // Only points around median +- cut filtered
1627 if (!graph) return 0;
1629 Int_t npoints0 = graph->GetN();
1632 Double_t *outx=new Double_t[npoints0];
1633 Double_t *outy=new Double_t[npoints0];
1636 if (npoints0<kMinPoints) return 0;
1637 for (Int_t iter=0; iter<3; iter++){
1639 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1640 if (graph->GetY()[ipoint]==0) continue;
1641 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
1642 outx[npoints] = graph->GetX()[ipoint];
1643 outy[npoints] = graph->GetY()[ipoint];
1646 if (npoints<=1) break;
1647 medianY =TMath::Median(npoints,outy);
1648 rmsY =TMath::RMS(npoints,outy);
1651 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1657 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
1659 // filter outlyer measurement
1660 // Only points with normalized errors median +- sigmaCut filtered
1662 Int_t kMinPoints=10;
1663 Int_t npoints0 = graph->GetN();
1665 Float_t medianErr=0, rmsErr=0;
1666 Double_t *outx=new Double_t[npoints0];
1667 Double_t *outy=new Double_t[npoints0];
1668 Double_t *erry=new Double_t[npoints0];
1669 Double_t *nerry=new Double_t[npoints0];
1670 Double_t *errx=new Double_t[npoints0];
1673 if (npoints0<kMinPoints) return 0;
1674 for (Int_t iter=0; iter<3; iter++){
1676 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1677 nerry[npoints] = graph->GetErrorY(ipoint);
1678 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
1679 erry[npoints] = graph->GetErrorY(ipoint);
1680 outx[npoints] = graph->GetX()[ipoint];
1681 outy[npoints] = graph->GetY()[ipoint];
1682 errx[npoints] = graph->GetErrorY(ipoint);
1685 if (npoints==0) break;
1686 medianErr=TMath::Median(npoints,erry);
1687 medianY =TMath::Median(npoints,outy);
1688 rmsErr =TMath::RMS(npoints,erry);
1690 TGraphErrors *graphOut=0;
1691 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
1700 void AliTPCcalibDButil::Sort(TGraph *graph){
1702 // sort array - neccessay for approx
1704 Int_t npoints = graph->GetN();
1705 Int_t *indexes=new Int_t[npoints];
1706 Double_t *outx=new Double_t[npoints];
1707 Double_t *outy=new Double_t[npoints];
1708 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
1709 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
1710 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
1711 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
1712 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
1715 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
1717 // smmoth graph - mean on the interval
1720 Int_t npoints = graph->GetN();
1721 Double_t *outy=new Double_t[npoints];
1723 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1724 Double_t lx=graph->GetX()[ipoint];
1725 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1726 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1727 if (index0<0) index0=0;
1728 if (index1>=npoints-1) index1=npoints-1;
1729 if ((index1-index0)>1){
1730 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1732 outy[ipoint]=graph->GetY()[ipoint];
1735 // TLinearFitter fitter(3,"pol2");
1736 // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1737 // Double_t lx=graph->GetX()[ipoint];
1738 // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1739 // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1740 // if (index0<0) index0=0;
1741 // if (index1>=npoints-1) index1=npoints-1;
1742 // fitter.ClearPoints();
1743 // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
1744 // if ((index1-index0)>1){
1745 // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1747 // outy[ipoint]=graph->GetY()[ipoint];
1753 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1754 graph->GetY()[ipoint] = outy[ipoint];
1759 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
1761 // Use constant interpolation outside of range
1764 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
1767 if (graph->GetN()<1){
1768 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
1771 if (xref<graph->GetX()[0]) return graph->GetY()[0];
1772 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
1773 return graph->Eval( xref);
1776 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
1778 // Filter DCS sensor information
1779 // ymin - minimal value
1781 // maxdy - maximal deirivative
1782 // sigmaCut - cut on values and derivative in terms of RMS distribution
1783 // Return value - accepted fraction
1787 // 0. Calculate median and rms of values in specified range
1788 // 1. Filter out outliers - median+-sigmaCut*rms
1789 // values replaced by median
1791 AliSplineFit * fit = sensor->GetFit();
1792 if (!fit) return 0.;
1793 Int_t nknots = fit->GetKnots();
1800 Double_t *yin0 = new Double_t[nknots];
1801 Double_t *yin1 = new Double_t[nknots];
1804 for (Int_t iknot=0; iknot< nknots; iknot++){
1805 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
1806 yin0[naccept] = fit->GetY0()[iknot];
1807 yin1[naccept] = fit->GetY1()[iknot];
1808 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
1818 Double_t medianY0=0, medianY1=0;
1819 Double_t rmsY0 =0, rmsY1=0;
1820 medianY0 = TMath::Median(naccept, yin0);
1821 medianY1 = TMath::Median(naccept, yin1);
1822 rmsY0 = TMath::RMS(naccept, yin0);
1823 rmsY1 = TMath::RMS(naccept, yin1);
1826 // 1. Filter out outliers - median+-sigmaCut*rms
1827 // values replaced by median
1828 // if replaced the derivative set to 0
1830 for (Int_t iknot=0; iknot< nknots; iknot++){
1832 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
1833 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
1834 if (nknots<2) fit->GetY1()[iknot]=0;
1835 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
1837 fit->GetY0()[iknot]=medianY0;
1838 fit->GetY1()[iknot]=0;
1845 return Float_t(naccept)/Float_t(nknots);
1848 Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
1850 // Filter temperature array
1851 // tempArray - array of temperatures -
1852 // ymin - minimal accepted temperature - default 15
1853 // ymax - maximal accepted temperature - default 22
1854 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
1855 // return value - fraction of filtered sensors
1856 const Double_t kMaxDy=0.1;
1857 Int_t nsensors=tempArray->NumSensors();
1858 if (nsensors==0) return 0.;
1860 for (Int_t isensor=0; isensor<nsensors; isensor++){
1861 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
1862 if (!sensor) continue;
1863 //printf("%d\n",isensor);
1864 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
1865 if (sensor->GetFit()==0){
1867 tempArray->RemoveSensorNum(isensor);
1872 return Float_t(naccept)/Float_t(nsensors);
1876 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
1879 // Input parameters:
1880 // deltaT - smoothing window (in seconds)
1881 // cutAbs - max distance of the time info to the median (in time bins)
1882 // cutSigma - max distance (in the RMS)
1883 // pcstream - optional debug streamer to store original and filtered info
1884 // Hardwired parameters:
1885 // kMinPoints =10; // minimal number of points to define the CE
1886 // kMinSectors=12; // minimal number of sectors to define sideCE
1888 // 0. Filter almost emty graphs (kMinPoints=10)
1889 // 1. calculate median and RMS per side
1890 // 2. Filter graphs - in respect with side medians
1891 // - cutAbs and cutDelta used
1892 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
1893 // 4. Calculate mean for A side and C side
1895 const Int_t kMinPoints =10; // minimal number of points to define the CE
1896 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
1897 const Int_t kMinTime =400; // minimal arrival time of CE
1898 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
1900 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1901 if (!cearray) return;
1906 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
1907 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernPressure");
1908 if ( tempMapCE && cavernPressureCE){
1910 Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
1911 FilterSensor(cavernPressureCE,960,1050,10, 5.);
1912 if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
1914 // recalculate P/T correction map for time of the CE
1915 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
1916 driftCalib->SetName("driftPTCE");
1917 driftCalib->SetTitle("driftPTCE");
1918 cearray->AddLast(driftCalib);
1922 // 0. Filter almost emty graphs
1925 for (Int_t i=0; i<72;i++){
1926 TGraph *graph= (TGraph*)arrT->At(i);
1927 if (!graph) continue;
1928 if (graph->GetN()<kMinPoints){
1930 delete graph; // delete empty graph
1933 if (tmin<0) tmin = graph->GetX()[0];
1934 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
1936 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
1937 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
1940 // 1. calculate median and RMS per side
1942 TArrayF arrA(100000), arrC(100000);
1944 Double_t medianA=0, medianC=0;
1945 Double_t rmsA=0, rmsC=0;
1946 for (Int_t isec=0; isec<72;isec++){
1947 TGraph *graph= (TGraph*)arrT->At(isec);
1948 if (!graph) continue;
1949 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1950 if (graph->GetY()[ipoint]<kMinTime) continue;
1951 if (nA>=arrA.fN) arrA.Set(nA*2);
1952 if (nC>=arrC.fN) arrC.Set(nC*2);
1953 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
1954 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
1958 medianA=TMath::Median(nA,arrA.fArray);
1959 rmsA =TMath::RMS(nA,arrA.fArray);
1962 medianC=TMath::Median(nC,arrC.fArray);
1963 rmsC =TMath::RMS(nC,arrC.fArray);
1966 // 2. Filter graphs - in respect with side medians
1968 TArrayD vecX(100000), vecY(100000);
1969 for (Int_t isec=0; isec<72;isec++){
1970 TGraph *graph= (TGraph*)arrT->At(isec);
1971 if (!graph) continue;
1972 Double_t median = (isec%36<18) ? medianA: medianC;
1973 Double_t rms = (isec%36<18) ? rmsA: rmsC;
1975 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1976 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
1977 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
1978 vecX[naccept]= graph->GetX()[ipoint];
1979 vecY[naccept]= graph->GetY()[ipoint];
1982 if (naccept<kMinPoints){
1983 arrT->AddAt(0,isec);
1984 delete graph; // delete empty graph
1987 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
1989 arrT->AddAt(graph2,isec);
1992 // 3. Cut in respect wit the graph median
1994 for (Int_t i=0; i<72;i++){
1995 TGraph *graph= (TGraph*)arrT->At(i);
1996 if (!graph) continue;
2000 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2001 if (!graphTS0) continue;
2002 if (graphTS0->GetN()<kMinPoints) {
2008 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2010 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2012 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2013 (*pcstream)<<"filterCE"<<
2018 "graphTS0.="<<graphTS0<<
2019 "graphTS.="<<graphTS<<
2023 if (!graphTS) continue;
2024 arrT->AddAt(graphTS,i);
2028 // Recalculate the mean time A side C side
2030 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2031 Int_t meanPoints=(nA+nC)/72; // mean number of points
2032 for (Int_t itime=0; itime<200; itime++){
2034 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2035 for (Int_t i=0; i<72;i++){
2036 TGraph *graph= (TGraph*)arrT->At(i);
2037 if (!graph) continue;
2038 if (graph->GetN()<(meanPoints/4)) continue;
2039 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2040 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2044 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2045 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2046 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2047 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2050 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2051 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2053 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2054 (*pcstream)<<"filterAC"<<
2063 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2064 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2065 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2066 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2067 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2068 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2071 if (nA<kMinSectors) arrT->AddAt(0,72);
2072 else arrT->AddAt(graphTSA,72);
2073 if (nC<kMinSectors) arrT->AddAt(0,73);
2074 else arrT->AddAt(graphTSC,73);
2078 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
2080 // Filter Drift velocity measurement using the tracks
2081 // 0. remove outlyers - error based
2085 const Int_t kMinPoints=1; // minimal number of points to define value
2086 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2089 for (Int_t i=0; i<arrT->GetEntries();i++){
2090 TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2091 if (!graph) continue;
2092 if (graph->GetN()<kMinPoints){
2097 TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2099 delete graph; arrT->AddAt(0,i); continue;
2101 if (graph2->GetN()<1) {
2102 delete graph; arrT->AddAt(0,i); continue;
2104 graph2->SetName(graph->GetName());
2105 graph2->SetTitle(graph->GetTitle());
2106 arrT->AddAt(graph2,i);
2108 (*pcstream)<<"filterTracks"<<
2113 "graph2.="<<graph2<<
2124 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2127 // get laser time offset
2128 // median around timeStamp+-deltaT
2129 // QA - chi2 needed for later usage - to be added
2130 // - currently cut on error
2133 Double_t kMinDelay=0.01;
2134 Double_t kMinDelayErr=0.0001;
2136 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2137 if (!array) return 0;
2138 TGraphErrors *tlaser=0;
2140 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2141 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2143 if (!tlaser) return 0;
2144 Int_t npoints0= tlaser->GetN();
2145 if (npoints0==0) return 0;
2146 Double_t *xlaser = new Double_t[npoints0];
2147 Double_t *ylaser = new Double_t[npoints0];
2149 for (Int_t i=0;i<npoints0;i++){
2151 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2152 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2153 xlaser[npoints]=tlaser->GetX()[npoints];
2154 ylaser[npoints]=tlaser->GetY()[npoints];
2159 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2160 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2161 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2162 if (index0<0) index0=0;
2163 if (index1>=npoints-1) index1=npoints-1;
2164 if (index1-index0<kMinPoints) return 0;
2166 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2167 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2176 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector *pcstream){
2178 // Filter Goofie data
2179 // goofieArray - points will be filtered
2180 // deltaT - smmothing time window
2181 // cutSigma - outler sigma cut in rms
2182 // minVn, maxVd- range absolute cut for variable vd/pt
2185 // Ignore goofie if not enough points
2187 const Int_t kMinPoints = 3;
2190 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2191 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2192 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2193 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2194 if (!graphvd) return;
2195 if (graphvd->GetN()<kMinPoints){
2197 goofieArray->GetSensorNum(2)->SetGraph(0);
2201 // 1. Caluclate medians of critical variables
2207 Double_t medianpt=0;
2208 Double_t medianvd=0, sigmavd=0;
2209 Double_t medianan=0;
2210 Double_t medianaf=0;
2211 Int_t entries=graphvd->GetN();
2212 Double_t yvdn[10000];
2215 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2216 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2217 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2218 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2219 yvdn[nvd++]=graphvd->GetY()[ipoint];
2221 if (nvd<kMinPoints){
2223 goofieArray->GetSensorNum(2)->SetGraph(0);
2227 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2228 if (nuni>=kMinPoints){
2229 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2231 medianvd = TMath::Median(nvd, yvdn);
2234 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2235 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2236 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2237 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2238 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2239 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2247 // 2. Make outlyer graph
2250 TGraph graphOut(*graphvd);
2251 for (Int_t i=0; i<entries;i++){
2253 Bool_t isOut=kFALSE;
2254 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2255 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2257 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2259 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2260 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2261 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2262 graphOut.GetY()[i]= (isOut)?1:0;
2265 if (nOK<kMinPoints) {
2267 goofieArray->GetSensorNum(2)->SetGraph(0);
2271 // 3. Filter out outlyers - and smooth
2273 TVectorF vmedianArray(goofieArray->NumSensors());
2274 TVectorF vrmsArray(goofieArray->NumSensors());
2275 Double_t xnew[10000];
2276 Double_t ynew[10000];
2278 junk.SetOwner(kTRUE);
2282 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2284 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2285 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2287 if (!sensor) continue;
2288 graphOld = sensor->GetGraph();
2290 sensor->SetGraph(0);
2292 for (Int_t i=0;i<entries;i++){
2293 if (graphOut.GetY()[i]>0.5) continue;
2294 xnew[nused]=graphOld->GetX()[i];
2295 ynew[nused]=graphOld->GetY()[i];
2298 graphNew = new TGraph(nused,xnew,ynew);
2299 junk.AddLast(graphNew);
2300 junk.AddLast(graphOld);
2302 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2304 junk.AddLast(graphNew0);
2305 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2307 junk.AddLast(graphNew1);
2308 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2310 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2311 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2312 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2313 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2314 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2315 vmedianArray[isensor]=median;
2321 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2322 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2323 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2324 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2325 (*pcstream)<<"goofieA"<<
2326 Form("isOK_%d.=",isensor)<<isOK<<
2327 Form("s_%d.=",isensor)<<sensor<<
2328 Form("gr_%d.=",isensor)<<graphOld<<
2329 Form("gr0_%d.=",isensor)<<graphNew0<<
2330 Form("gr1_%d.=",isensor)<<graphNew1<<
2331 Form("gr2_%d.=",isensor)<<graphNew2;
2332 if (isOK) sensor->SetGraph(graphNew2);
2334 (*pcstream)<<"goofieA"<<
2335 "vmed.="<<&vmedianArray<<
2336 "vrms.="<<&vrmsArray<<
2338 junk.Delete(); // delete temoprary graphs