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"
56 #include "AliRelAlignerKalman.h"
58 const Float_t kAlmost0=1.e-30;
60 ClassImp(AliTPCcalibDButil)
61 AliTPCcalibDButil::AliTPCcalibDButil() :
69 fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
81 fRefPedestalMasked(0x0),
85 fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
86 fRefPulserMasked(0x0),
93 fRefALTROAcqStart(0x0),
94 fRefALTROAcqStop(0x0),
99 fMapper(new AliTPCmapper(0x0)),
100 fNpulserOutliers(-1),
102 fCETmaxLimitAbs(1.5),
103 fPulTmaxLimitAbs(1.5),
106 fRuns(0), // run list with OCDB info
107 fRunsStart(0), // start time for given run
108 fRunsStop(0) // stop time for given run
114 //_____________________________________________________________________________________
115 AliTPCcalibDButil::~AliTPCcalibDButil()
120 delete fPulserOutlier;
121 delete fRefPulserOutlier;
123 if (fRefPadNoise) delete fRefPadNoise;
124 if (fRefPedestals) delete fRefPedestals;
125 if (fRefPedestalMasked) delete fRefPedestalMasked;
126 if (fRefPulserTmean) delete fRefPulserTmean;
127 if (fRefPulserTrms) delete fRefPulserTrms;
128 if (fRefPulserQmean) delete fRefPulserQmean;
129 if (fRefPulserMasked) delete fRefPulserMasked;
130 if (fRefCETmean) delete fRefCETmean;
131 if (fRefCETrms) delete fRefCETrms;
132 if (fRefCEQmean) delete fRefCEQmean;
133 if (fRefCEMasked) delete fRefCEMasked;
134 if (fRefALTROFPED) delete fRefALTROFPED;
135 if (fRefALTROZsThr) delete fRefALTROZsThr;
136 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
137 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
138 if (fRefALTROMasked) delete fRefALTROMasked;
139 if (fRefCalibRaw) delete fRefCalibRaw;
140 if (fCurrentRefMap) delete fCurrentRefMap;
142 //_____________________________________________________________________________________
143 void AliTPCcalibDButil::UpdateFromCalibDB()
146 // Update pointers from calibDB
148 if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
149 fCalibDB->UpdateNonRec(); // load all infromation now
150 fPadNoise=fCalibDB->GetPadNoise();
151 fPedestals=fCalibDB->GetPedestals();
152 fPulserTmean=fCalibDB->GetPulserTmean();
153 fPulserTrms=fCalibDB->GetPulserTrms();
154 fPulserQmean=fCalibDB->GetPulserQmean();
155 fCETmean=fCalibDB->GetCETmean();
156 fCETrms=fCalibDB->GetCETrms();
157 fCEQmean=fCalibDB->GetCEQmean();
158 fALTROMasked=fCalibDB->GetALTROMasked();
159 fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
160 fCalibRaw=fCalibDB->GetCalibRaw();
161 fDataQA=fCalibDB->GetDataQA();
162 UpdatePulserOutlierMap();
163 // SetReferenceRun();
164 // UpdateRefDataFromOCDB();
166 //_____________________________________________________________________________________
167 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
168 Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
171 // Process the CE data for this run
172 // the return TVectorD arrays contian the results of the fit
173 // noutliersCE contains the number of pads marked as outliers,
174 // not including masked and edge pads
177 //retrieve CE and ALTRO data
179 TString fitString(fitFormula);
180 fitString.ReplaceAll("++","#");
181 Int_t ndim=fitString.CountChar('#')+2;
182 fitResultsA.ResizeTo(ndim);
183 fitResultsC.ResizeTo(ndim);
192 if (outCE) out=outCE;
193 else out=new AliTPCCalPad("outCE","outCE");
194 AliTPCCalROC *rocMasked=0x0;
195 //loop over all channels
196 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
197 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
198 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
199 AliTPCCalROC *rocOut=out->GetCalROC(iroc);
201 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
205 //add time offset to IROCs
206 if (iroc<AliTPCROC::Instance()->GetNInnerSector())
207 rocData->Add(fIrocTimeOffset);
209 UInt_t nrows=rocData->GetNrows();
210 for (UInt_t irow=0;irow<nrows;++irow){
211 UInt_t npads=rocData->GetNPads(irow);
212 for (UInt_t ipad=0;ipad<npads;++ipad){
213 rocOut->SetValue(irow,ipad,0);
214 //exclude masked pads
215 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
216 rocOut->SetValue(irow,ipad,1);
219 //exclude first two rows in IROC and last two rows in OROC
221 if (irow<2) rocOut->SetValue(irow,ipad,1);
223 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
226 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
227 Float_t valTmean=rocData->GetValue(irow,ipad);
228 //exclude values that are exactly 0
229 if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
230 rocOut->SetValue(irow,ipad,1);
233 // exclude channels with too large variations
234 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
235 rocOut->SetValue(irow,ipad,1);
243 Float_t chi2Af,chi2Cf;
244 fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
247 if (!outCE) delete out;
249 //_____________________________________________________________________________________
250 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
251 TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
252 Float_t &driftTimeA, Float_t &driftTimeC )
255 // Calculate statistical information from the CE graphs for drift time and charge
259 vecTEntries.ResizeTo(72);
260 vecTMean.ResizeTo(72);
261 vecTRMS.ResizeTo(72);
262 vecTMedian.ResizeTo(72);
263 vecQEntries.ResizeTo(72);
264 vecQMean.ResizeTo(72);
265 vecQRMS.ResizeTo(72);
266 vecQMedian.ResizeTo(72);
277 TObjArray *arrT=fCalibDB->GetCErocTtime();
278 TObjArray *arrQ=fCalibDB->GetCErocQtime();
280 for (Int_t isec=0;isec<74;++isec){
281 TGraph *gr=(TGraph*)arrT->At(isec);
284 Int_t npoints = gr->GetN();
285 values.ResizeTo(npoints);
287 //skip first points, theres always some problems with finding the CE position
288 for (Int_t ipoint=4; ipoint<npoints; ipoint++){
289 if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
290 values[nused]=gr->GetY()[ipoint];
295 if (isec<72) vecTEntries[isec]= nused;
298 vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
299 vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
300 vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
301 } else if (isec==72){
302 driftTimeA=TMath::Median(nused,values.GetMatrixArray());
303 } else if (isec==73){
304 driftTimeC=TMath::Median(nused,values.GetMatrixArray());
310 for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
311 TGraph *gr=(TGraph*)arrQ->At(isec);
314 Int_t npoints = gr->GetN();
315 values.ResizeTo(npoints);
317 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
318 if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
319 values[nused]=gr->GetY()[ipoint];
324 vecQEntries[isec]= nused;
326 vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
327 vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray());
328 vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray());
334 //_____________________________________________________________________________________
335 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
336 TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
337 Int_t &nonMaskedZero, Int_t &nNaN)
340 // process noise data
341 // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
342 // OROCs small pads [2] and OROCs large pads [3]
343 // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
344 // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
347 //set proper size and reset
348 const UInt_t infoSize=4;
349 vNoiseMean.ResizeTo(infoSize);
350 vNoiseMeanSenRegions.ResizeTo(infoSize);
351 vNoiseRMS.ResizeTo(infoSize);
352 vNoiseRMSSenRegions.ResizeTo(infoSize);
354 vNoiseMeanSenRegions.Zero();
356 vNoiseRMSSenRegions.Zero();
360 TVectorD c(infoSize);
361 TVectorD cs(infoSize);
365 //retrieve noise and ALTRO data
366 if (!fPadNoise) return;
367 AliTPCCalROC *rocMasked=0x0;
368 //create IROC, OROC1, OROC2 and sensitive region masks
369 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
370 AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
371 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
372 UInt_t nrows=noiseROC->GetNrows();
373 for (UInt_t irow=0;irow<nrows;++irow){
374 UInt_t npads=noiseROC->GetNPads(irow);
375 for (UInt_t ipad=0;ipad<npads;++ipad){
376 //don't use masked channels;
377 if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
378 Float_t noiseVal=noiseROC->GetValue(irow,ipad);
380 if (noiseVal<kAlmost0) {
385 if ( !(noiseVal<10000000) ){
386 // printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
390 Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
391 Int_t masksen=1; // sensitive pards are not masked (0)
392 if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
393 if (isec<AliTPCROC::Instance()->GetNInnerSector()){
395 if (irow>19&&irow<46){
396 if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
399 vNoiseMean[type]+=noiseVal;
400 vNoiseRMS[type]+=noiseVal*noiseVal;
403 vNoiseMeanSenRegions[type]+=noiseVal;
404 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
409 //define sensive regions
410 if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
412 Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
413 if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
415 if ((Int_t)irow<par.GetNRowUp1()){
418 vNoiseMean[type]+=noiseVal;
419 vNoiseRMS[type]+=noiseVal*noiseVal;
422 vNoiseMeanSenRegions[type]+=noiseVal;
423 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
429 vNoiseMean[type]+=noiseVal;
430 vNoiseRMS[type]+=noiseVal*noiseVal;
433 vNoiseMeanSenRegions[type]+=noiseVal;
434 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
441 vNoiseMean[type]+=noiseVal;
442 vNoiseRMS[type]+=noiseVal*noiseVal;
445 vNoiseMeanSenRegions[type]+=noiseVal;
446 vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
451 }//end loop sectors (rocs)
453 //calculate mean and RMS
454 const Double_t verySmall=0.0000000001;
455 for (UInt_t i=0;i<infoSize;++i){
462 // printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
463 mean=vNoiseMean[i]/c[i];
465 rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
470 if (cs[i]>verySmall){
471 meanSen=vNoiseMeanSenRegions[i]/cs[i];
472 rmsSen=vNoiseRMSSenRegions[i];
473 rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
475 vNoiseMeanSenRegions[i]=meanSen;
476 vNoiseRMSSenRegions[i]=rmsSen;
480 //_____________________________________________________________________________________
481 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
484 // Process the Pulser information
485 // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
488 const UInt_t infoSize=4;
489 //reset counters to error number
490 vMeanTime.ResizeTo(infoSize);
493 TVectorD c(infoSize);
494 //retrieve pulser and ALTRO data
495 if (!fPulserTmean) return;
498 AliTPCCalROC *rocOut=0x0;
499 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
500 AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
501 if (!tmeanROC) continue;
502 rocOut=fPulserOutlier->GetCalROC(isec);
503 UInt_t nchannels=tmeanROC->GetNchannels();
504 for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
505 if (rocOut && rocOut->GetValue(ichannel)) continue;
506 Float_t val=tmeanROC->GetValue(ichannel);
508 vMeanTime[type]+=val;
513 for (UInt_t itype=0; itype<infoSize; ++itype){
514 if (c[itype]>0) vMeanTime[itype]/=c[itype];
515 else vMeanTime[itype]=0;
518 //_____________________________________________________________________________________
519 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
522 // Get Values from ALTRO configuration data
525 if (!fALTROMasked) return;
527 for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
528 AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
529 for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
530 if (rocMasked->GetValue(ichannel)) ++nMasked;
534 //_____________________________________________________________________________________
535 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
538 // Proces Goofie values, return statistical information of the currently set goofieArray
539 // The meaning of the entries are given below
541 1 TPC_ANODE_I_A00_STAT
543 3 TPC_DVM_DriftVelocity
548 8 TPC_DVM_NumberOfSparks
549 9 TPC_DVM_PeakAreaFar
550 10 TPC_DVM_PeakAreaNear
551 11 TPC_DVM_PeakPosFar
552 12 TPC_DVM_PeakPosNear
558 18 TPC_DVM_TemperatureS1
562 vecEntries.ResizeTo(nsensors);
563 vecMedian.ResizeTo(nsensors);
564 vecMean.ResizeTo(nsensors);
565 vecRMS.ResizeTo(nsensors);
572 Double_t kEpsilon=0.0000000001;
573 Double_t kBig=100000000000.;
574 Int_t nsensors = fGoofieArray->NumSensors();
575 vecEntries.ResizeTo(nsensors);
576 vecMedian.ResizeTo(nsensors);
577 vecMean.ResizeTo(nsensors);
578 vecRMS.ResizeTo(nsensors);
580 for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
581 AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
582 if (gsensor && gsensor->GetGraph()){
583 Int_t npoints = gsensor->GetGraph()->GetN();
585 values.ResizeTo(npoints);
587 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
588 if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
589 TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
590 values[nused]=gsensor->GetGraph()->GetY()[ipoint];
595 vecEntries[isensor]= nused;
597 vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
598 vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray());
599 vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray());
604 //_____________________________________________________________________________________
605 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
608 // check the variations of the pedestal data to the reference pedestal data
609 // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
612 TVectorF vThres(npar); //thresholds
613 Int_t nActive=0; //number of active channels
615 //reset and set thresholds
616 pedestalDeviations.ResizeTo(npar);
617 for (Int_t i=0;i<npar;++i){
618 pedestalDeviations.GetMatrixArray()[i]=0;
619 vThres.GetMatrixArray()[i]=(i+1)*.5;
621 //check all needed data is available
622 if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
623 //loop over all channels
624 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
625 AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
626 AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
627 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
628 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
629 UInt_t nrows=mROC->GetNrows();
630 for (UInt_t irow=0;irow<nrows;++irow){
631 UInt_t npads=mROC->GetNPads(irow);
632 for (UInt_t ipad=0;ipad<npads;++ipad){
633 //don't use masked channels;
634 if (mROC ->GetValue(irow,ipad)) continue;
635 if (mRefROC->GetValue(irow,ipad)) continue;
636 Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
637 for (Int_t i=0;i<npar;++i){
638 if (deviation>vThres[i])
639 ++pedestalDeviations.GetMatrixArray()[i];
646 for (Int_t i=0;i<npar;++i){
647 pedestalDeviations.GetMatrixArray()[i]/=nActive;
651 //_____________________________________________________________________________________
652 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
655 // check the variations of the noise data to the reference noise data
656 // thresholds are 5, 10, 15 and 20 percent respectively.
659 TVectorF vThres(npar); //thresholds
660 Int_t nActive=0; //number of active channels
662 //reset and set thresholds
663 noiseDeviations.ResizeTo(npar);
664 for (Int_t i=0;i<npar;++i){
665 noiseDeviations.GetMatrixArray()[i]=0;
666 vThres.GetMatrixArray()[i]=(i+1)*.05;
668 //check all needed data is available
669 if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
670 //loop over all channels
671 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
672 AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
673 AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
674 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
675 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
676 UInt_t nrows=mROC->GetNrows();
677 for (UInt_t irow=0;irow<nrows;++irow){
678 UInt_t npads=mROC->GetNPads(irow);
679 for (UInt_t ipad=0;ipad<npads;++ipad){
680 //don't use masked channels;
681 if (mROC ->GetValue(irow,ipad)) continue;
682 if (mRefROC->GetValue(irow,ipad)) continue;
683 if (nRefROC->GetValue(irow,ipad)==0) continue;
684 Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
685 for (Int_t i=0;i<npar;++i){
686 if (deviation>vThres[i])
687 ++noiseDeviations.GetMatrixArray()[i];
694 for (Int_t i=0;i<npar;++i){
695 noiseDeviations.GetMatrixArray()[i]/=nActive;
699 //_____________________________________________________________________________________
700 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
701 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
704 // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
705 // thresholds are .5, 1, 5 and 10 percent respectively.
709 TVectorF vThres(npar); //thresholds
710 Int_t nActive=0; //number of active channels
712 //reset and set thresholds
713 pulserQdeviations.ResizeTo(npar);
714 for (Int_t i=0;i<npar;++i){
715 pulserQdeviations.GetMatrixArray()[i]=0;
720 vThres.GetMatrixArray()[0]=.005;
721 vThres.GetMatrixArray()[1]=.01;
722 vThres.GetMatrixArray()[2]=.05;
723 vThres.GetMatrixArray()[3]=.1;
724 //check all needed data is available
725 if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
727 UpdateRefPulserOutlierMap();
728 //loop over all channels
729 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
730 AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
731 AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
732 AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
733 // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
734 AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
735 AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
736 AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
737 Float_t ptmean=ptROC->GetMean(oROC);
738 UInt_t nrows=mROC->GetNrows();
739 for (UInt_t irow=0;irow<nrows;++irow){
740 UInt_t npads=mROC->GetNPads(irow);
741 for (UInt_t ipad=0;ipad<npads;++ipad){
742 //don't use masked channels;
743 if (mROC ->GetValue(irow,ipad)) continue;
744 if (mRefROC->GetValue(irow,ipad)) continue;
745 //don't user edge pads
746 if (ipad==0||ipad==npads-1) continue;
748 Float_t pq=pqROC->GetValue(irow,ipad);
749 Float_t pqRef=pqRefROC->GetValue(irow,ipad);
750 Float_t pt=ptROC->GetValue(irow,ipad);
751 // Float_t ptRef=ptRefROC->GetValue(irow,ipad);
753 Float_t deviation=TMath::Abs(pq/pqRef-1);
754 for (Int_t i=0;i<npar;++i){
755 if (deviation>vThres[i])
756 ++pulserQdeviations.GetMatrixArray()[i];
758 if (pqRef>11&&pq<11) ++npadsOffAdd;
761 if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
767 for (Int_t i=0;i<npar;++i){
768 pulserQdeviations.GetMatrixArray()[i]/=nActive;
773 //_____________________________________________________________________________________
774 void AliTPCcalibDButil::UpdatePulserOutlierMap()
777 // Update the outlier map of the pulser data
779 PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
781 //_____________________________________________________________________________________
782 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
785 // Update the outlier map of the pulser reference data
787 PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
789 //_____________________________________________________________________________________
790 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
793 // Create a map that contains outliers from the Pulser calibration data.
794 // The outliers include masked channels, edge pads and pads with
795 // too large timing and charge variations.
796 // fNpulserOutliers is the number of outliers in the Pulser calibration data.
797 // those do not contain masked and edge pads
801 pulOut->Multiply(0.);
805 AliTPCCalROC *rocMasked=0x0;
809 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
810 AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
811 AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
812 AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
813 if (!tmeanROC||!qmeanROC) {
814 //reset outliers in this ROC
815 outROC->Multiply(0.);
818 if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
820 // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
821 // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
822 UInt_t nrows=tmeanROC->GetNrows();
823 for (UInt_t irow=0;irow<nrows;++irow){
824 UInt_t npads=tmeanROC->GetNPads(irow);
825 for (UInt_t ipad=0;ipad<npads;++ipad){
826 Int_t outlier=0,masked=0;
827 Float_t q=qmeanROC->GetValue(irow,ipad);
828 Float_t t=tmeanROC->GetValue(irow,ipad);
829 //masked channels are outliers
830 if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
831 //edge pads are outliers
832 if (ipad==0||ipad==npads-1) masked=1;
833 //channels with too large charge or timing deviation from the meadian are outliers
834 // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
835 if (q<fPulQminLimit && !masked) outlier=1;
837 if ( !(q<10000000) || !(t<10000000)) outlier=1;
838 outROC->SetValue(irow,ipad,outlier+masked);
839 fNpulserOutliers+=outlier;
844 //_____________________________________________________________________________________
845 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
848 // Create pad time0 object from pulser and/or CE data, depending on the selected model
849 // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
850 // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
851 // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
853 // In case model 2 is invoked - gy arival time gradient is also returned
857 AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
858 // decide between different models
859 if (model==0||model==1){
861 if (model==1) ProcessPulser(vMean);
862 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
863 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
864 if (!rocPulTmean) continue;
865 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
866 AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
867 Float_t mean=rocPulTmean->GetMean(rocOut);
868 //treat case where a whole partition is masked
869 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
874 UInt_t nrows=rocTime0->GetNrows();
875 for (UInt_t irow=0;irow<nrows;++irow){
876 UInt_t npads=rocTime0->GetNPads(irow);
877 for (UInt_t ipad=0;ipad<npads;++ipad){
878 Float_t time=rocPulTmean->GetValue(irow,ipad);
879 //in case of an outlier pad use the mean of the altro values.
880 //This should be the most precise guess in that case.
881 if (rocOut->GetValue(irow,ipad)) {
882 time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
883 if ( TMath::Abs(time)<kAlmost0 ) time=mean;
885 Float_t val=time-mean;
886 rocTime0->SetValue(irow,ipad,val);
890 } else if (model==2){
891 Double_t pgya,pgyc,pchi2a,pchi2c;
892 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
893 fCETmean->Add(padPulser,-1.);
895 AliTPCCalPad outCE("outCE","outCE");
897 ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
898 AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
899 // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
900 if (!padFit) { delete padPulser; return 0;}
903 fCETmean->Add(padPulser,1.);
904 padTime0->Add(fCETmean);
905 padTime0->Add(padFit,-1);
910 for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
911 AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
912 AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
913 AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
914 AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
915 rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
916 AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
917 Float_t mean=rocPulTmean->GetMean(rocOutPul);
918 if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
919 UInt_t nrows=rocTime0->GetNrows();
920 for (UInt_t irow=0;irow<nrows;++irow){
921 UInt_t npads=rocTime0->GetNPads(irow);
922 for (UInt_t ipad=0;ipad<npads;++ipad){
923 Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
924 if (rocOutCE->GetValue(irow,ipad)){
925 Float_t valOut=rocCEfit->GetValue(irow,ipad);
926 if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
927 rocTime0->SetValue(irow,ipad,valOut);
935 Double_t median = padTime0->GetMedian();
936 padTime0->Add(-median); // normalize to median
939 //_____________________________________________________________________________________
940 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *const rocOut)
943 // GetMeanAlto information
945 if (roc==0) return 0.;
946 const Int_t sector=roc->GetSector();
947 AliTPCROC *tpcRoc=AliTPCROC::Instance();
948 const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
952 //loop over a small range around the requested pad (+-10 rows/pads)
953 for (Int_t irow=row-10;irow<row+10;++irow){
954 if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
955 for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
956 if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
957 const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
958 if (altroRoc!=altroCurr) continue;
959 if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
960 Float_t val=roc->GetValue(irow,ipad);
968 //_____________________________________________________________________________________
969 void AliTPCcalibDButil::SetRefFile(const char* filename)
972 // load cal pad objects form the reference file
974 TDirectory *currDir=gDirectory;
976 fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
977 fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
979 fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
980 fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
981 fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
983 fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
984 fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
985 fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
987 // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
988 // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
989 // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
990 // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
991 fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
995 //_____________________________________________________________________________________
996 void AliTPCcalibDButil::UpdateRefDataFromOCDB()
999 // set reference data from OCDB Reference map
1002 AliWarning("Referenc map not set!");
1007 AliCDBEntry* entry = 0x0;
1008 Bool_t hasAnyChanged=kFALSE;
1011 cdbPath="TPC/Calib/Pedestals";
1012 if (HasRefChanged(cdbPath.Data())){
1013 hasAnyChanged=kTRUE;
1014 //delete old entries
1015 if (fRefPedestals) delete fRefPedestals;
1016 if (fRefPedestalMasked) delete fRefPedestalMasked;
1017 fRefPedestals=fRefPedestalMasked=0x0;
1019 entry=GetRefEntry(cdbPath.Data());
1021 entry->SetOwner(kTRUE);
1022 fRefPedestals=GetRefCalPad(entry);
1024 fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1029 cdbPath="TPC/Calib/PadNoise";
1030 if (HasRefChanged(cdbPath.Data())){
1031 hasAnyChanged=kTRUE;
1033 if (fRefPadNoise) delete fRefPadNoise;
1036 entry=GetRefEntry(cdbPath.Data());
1038 entry->SetOwner(kTRUE);
1039 fRefPadNoise=GetRefCalPad(entry);
1045 cdbPath="TPC/Calib/Pulser";
1046 if (HasRefChanged(cdbPath.Data())){
1047 hasAnyChanged=kTRUE;
1048 //delete old entries
1049 if (fRefPulserTmean) delete fRefPulserTmean;
1050 if (fRefPulserTrms) delete fRefPulserTrms;
1051 if (fRefPulserQmean) delete fRefPulserQmean;
1052 if (fRefPulserMasked) delete fRefPulserMasked;
1053 fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1055 entry=GetRefEntry(cdbPath.Data());
1057 entry->SetOwner(kTRUE);
1058 fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1059 fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1060 fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1062 fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1067 cdbPath="TPC/Calib/CE";
1068 if (HasRefChanged(cdbPath.Data())){
1069 hasAnyChanged=kTRUE;
1070 //delete old entries
1071 if (fRefCETmean) delete fRefCETmean;
1072 if (fRefCETrms) delete fRefCETrms;
1073 if (fRefCEQmean) delete fRefCEQmean;
1074 if (fRefCEMasked) delete fRefCEMasked;
1075 fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1077 entry=GetRefEntry(cdbPath.Data());
1079 entry->SetOwner(kTRUE);
1080 fRefCETmean=GetRefCalPad(entry,"CETmean");
1081 fRefCETrms=GetRefCalPad(entry,"CETrms");
1082 fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1084 fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1089 cdbPath="TPC/Calib/AltroConfig";
1090 if (HasRefChanged(cdbPath.Data())){
1091 hasAnyChanged=kTRUE;
1092 //delete old entries
1093 if (fRefALTROFPED) delete fRefALTROFPED;
1094 if (fRefALTROZsThr) delete fRefALTROZsThr;
1095 if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1096 if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1097 if (fRefALTROMasked) delete fRefALTROMasked;
1098 fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1100 entry=GetRefEntry(cdbPath.Data());
1102 entry->SetOwner(kTRUE);
1103 fRefALTROFPED=GetRefCalPad(entry,"FPED");
1104 fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1105 fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1106 fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1107 fRefALTROMasked=GetRefCalPad(entry,"Masked");
1114 cdbPath="TPC/Calib/Raw";
1115 if (HasRefChanged(cdbPath.Data())){
1116 hasAnyChanged=kTRUE;
1118 if (fRefCalibRaw) delete fRefCalibRaw;
1120 entry=GetRefEntry(cdbPath.Data());
1122 entry->SetOwner(kTRUE);
1123 TObjArray *arr=(TObjArray*)entry->GetObject();
1125 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1127 fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1134 cdbPath="TPC/Calib/QA";
1135 if (HasRefChanged(cdbPath.Data())){
1136 hasAnyChanged=kTRUE;
1138 if (fRefDataQA) delete fRefDataQA;
1140 entry=GetRefEntry(cdbPath.Data());
1142 entry->SetOwner(kTRUE);
1143 fDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1145 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1147 fRefDataQA=(AliTPCdataQA*)fDataQA->Clone();
1154 //update current reference maps
1156 if (fCurrentRefMap) delete fCurrentRefMap;
1157 fCurrentRefMap=(TMap*)fRefMap->Clone();
1160 //_____________________________________________________________________________________
1161 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1164 // TObjArray object type case
1165 // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1167 AliTPCCalPad *pad=0x0;
1168 TObjArray *arr=(TObjArray*)entry->GetObject();
1170 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1173 pad=(AliTPCCalPad*)arr->FindObject(objName);
1175 AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
1178 return (AliTPCCalPad*)pad->Clone();
1180 //_____________________________________________________________________________________
1181 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1184 // AliTPCCalPad object type case
1185 // cast object to a calPad and store it in 'pad'
1187 AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1189 AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1192 pad=(AliTPCCalPad*)pad->Clone();
1195 //_____________________________________________________________________________________
1196 AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1199 // set altro masked channel map for 'cdbPath'
1201 AliTPCCalPad* pad=0x0;
1202 const Int_t run=GetReferenceRun(cdbPath);
1204 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1207 AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1209 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1212 pad=GetRefCalPad(entry,"Masked");
1213 if (pad) pad->SetNameTitle(name,name);
1214 entry->SetOwner(kTRUE);
1218 //_____________________________________________________________________________________
1219 void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1221 // Get Reference map
1223 if (run<0) run=fCalibDB->GetRun();
1224 TString cdbPath="TPC/Calib/Ref";
1225 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1227 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
1231 entry->SetOwner(kTRUE);
1232 fRefMap=(TMap*)(entry->GetObject());
1233 AliCDBId &id=entry->GetId();
1234 fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1236 //_____________________________________________________________________________________
1237 Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1240 // check whether a reference cdb entry has changed
1242 if (!fCurrentRefMap) return kTRUE;
1243 if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1246 //_____________________________________________________________________________________
1247 AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1250 // get the reference AliCDBEntry for 'cdbPath'
1252 const Int_t run=GetReferenceRun(cdbPath);
1254 AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1257 AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1259 AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1264 //_____________________________________________________________________________________
1265 Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
1267 // Get reference run number for the specified OCDB path
1269 if (!fCurrentRefMap) return -2;
1270 TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1271 if (!str) return -2;
1272 return (const Int_t)str->GetString().Atoi();
1274 //_____________________________________________________________________________________
1275 Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1277 // Get reference run number for the specified OCDB path
1279 if (!fRefMap) return -1;
1280 TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1281 if (!str) return -1;
1282 return (const Int_t)str->GetString().Atoi();
1284 //_____________________________________________________________________________________
1285 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){
1287 // Author: marian.ivanov@cern.ch
1289 // Create outlier map for CE study
1291 // Return value - outlyer map
1292 // noutlyersCE - number of outlyers
1293 // minSignal - minimal total Q signal
1294 // cutRMSMin - minimal width of the signal in respect to the median
1295 // cutRMSMax - maximal width of the signal in respect to the median
1296 // cutMaxDistT - maximal deviation from time median per chamber
1298 // Outlyers criteria:
1299 // 0. Exclude masked pads
1300 // 1. Exclude first two rows in IROC and last two rows in OROC
1301 // 2. Exclude edge pads
1302 // 3. Exclude channels with too large variations
1303 // 4. Exclude pads with too small signal
1304 // 5. Exclude signal with outlyers RMS
1305 // 6. Exclude channels to far from the chamber median
1307 //create outlier map
1308 AliTPCCalPad *out=ceOut;
1309 if (!out) out= new AliTPCCalPad("outCE","outCE");
1310 AliTPCCalROC *rocMasked=0x0;
1311 if (!fCETmean) return 0;
1312 if (!fCETrms) return 0;
1313 if (!fCEQmean) return 0;
1315 //loop over all channels
1317 Double_t rmsMedian = fCETrms->GetMedian();
1318 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1319 AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1320 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1321 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1322 AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc);
1323 AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc);
1324 Double_t trocMedian = rocData->GetMedian();
1327 noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1333 UInt_t nrows=rocData->GetNrows();
1334 for (UInt_t irow=0;irow<nrows;++irow){
1335 UInt_t npads=rocData->GetNPads(irow);
1336 for (UInt_t ipad=0;ipad<npads;++ipad){
1337 rocOut->SetValue(irow,ipad,0);
1338 Float_t valTmean=rocData->GetValue(irow,ipad);
1339 Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1340 Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1341 //0. exclude masked pads
1342 if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1343 rocOut->SetValue(irow,ipad,1);
1346 //1. exclude first two rows in IROC and last two rows in OROC
1348 if (irow<2) rocOut->SetValue(irow,ipad,1);
1350 if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1352 //2. exclude edge pads
1353 if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1354 //exclude values that are exactly 0
1355 if ( TMath::Abs(valTmean)<kAlmost0) {
1356 rocOut->SetValue(irow,ipad,1);
1359 //3. exclude channels with too large variations
1360 if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1361 rocOut->SetValue(irow,ipad,1);
1365 //4. exclude channels with too small signal
1366 if (valQmean<minSignal) {
1367 rocOut->SetValue(irow,ipad,1);
1371 //5. exclude channels with too small rms
1372 if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1373 rocOut->SetValue(irow,ipad,1);
1377 //6. exclude channels to far from the chamber median
1378 if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1379 rocOut->SetValue(irow,ipad,1);
1390 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1392 // Author: marian.ivanov@cern.ch
1394 // Create outlier map for Pulser
1396 // Return value - outlyer map
1397 // noutlyersPulser - number of outlyers
1398 // cutTime - absolute cut - distance to the median of chamber
1399 // cutnRMSQ - nsigma cut from median q distribution per chamber
1400 // cutnRMSrms - nsigma cut from median rms distribution
1401 // Outlyers criteria:
1402 // 0. Exclude masked pads
1403 // 1. Exclude time outlyers (default 3 time bins)
1404 // 2. Exclude q outlyers (default 5 sigma)
1405 // 3. Exclude rms outlyers (default 5 sigma)
1407 AliTPCCalPad *out=pulserOut;
1408 if (!out) out= new AliTPCCalPad("outPulser","outPulser");
1409 AliTPCCalROC *rocMasked=0x0;
1410 if (!fPulserTmean) return 0;
1411 if (!fPulserTrms) return 0;
1412 if (!fPulserQmean) return 0;
1414 //loop over all channels
1416 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1417 if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1418 AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc);
1419 AliTPCCalROC *rocOut = out->GetCalROC(iroc);
1420 AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc);
1421 AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1423 Double_t rocMedianT = rocData->GetMedian();
1424 Double_t rocMedianQ = rocPulserQ->GetMedian();
1425 Double_t rocRMSQ = rocPulserQ->GetRMS();
1426 Double_t rocMedianTrms = rocPulserTrms->GetMedian();
1427 Double_t rocRMSTrms = rocPulserTrms->GetRMS();
1428 for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1429 rocOut->SetValue(ichannel,0);
1430 Float_t valTmean=rocData->GetValue(ichannel);
1431 Float_t valQmean=rocPulserQ->GetValue(ichannel);
1432 Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1434 if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1435 if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1436 if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1437 rocOut->SetValue(ichannel,isOut);
1438 if (isOut) noutliersPulser++;
1445 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1447 // Author : Marian Ivanov
1448 // Create pad time0 correction map using information from the CE and from pulser
1451 // Return PadTime0 to be used for time0 relative alignment
1452 // if dump file specified intermediat results are dumped to the fiel and can be visualized
1453 // using $ALICE_ROOT/TPC/script/gui application
1455 // fitResultsA - fitParameters A side
1456 // fitResultsC - fitParameters C side
1457 // chi2A - chi2/ndf for A side (assuming error 1 time bin)
1458 // chi2C - chi2/ndf for C side (assuming error 1 time bin)
1462 // 1. Find outlier map for CE
1463 // 2. Find outlier map for Pulser
1464 // 3. Replace outlier by median at given sector (median without outliers)
1465 // 4. Substract from the CE data pulser
1466 // 5. Fit the CE with formula
1467 // 5.1) (IROC-OROC) offset
1471 // 5.5) (IROC-OROC)*(lx-xmid)
1473 // 6. Substract gy fit dependence from the CE data
1474 // 7. Add pulser back to CE data
1475 // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1476 // 9. return CE data
1478 // Time0 <= padCE = padCEin -padCEfitGy - if not outlier
1479 // Time0 <= padCE = padFitAll-padCEfitGy - if outlier
1482 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)";
1483 // output for fit formula
1484 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)";
1485 // gy part of formula
1486 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)";
1489 if (!fCETmean) return 0;
1490 Double_t pgya,pgyc,pchi2a,pchi2c;
1491 AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1492 AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut);
1494 AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1495 AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean);
1496 AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean);
1497 AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut");
1498 padPulser->SetName("padPulser");
1499 padPulserOut->SetName("padPulserOut");
1500 padCE->SetName("padCE");
1501 padCEIn->SetName("padCEIn");
1502 padCEOut->SetName("padCEOut");
1503 padOut->SetName("padOut");
1506 // make combined outlyers map
1507 // and replace outlyers in maps with median for chamber
1509 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1510 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1511 AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc);
1512 AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1513 AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc);
1514 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1515 Double_t ceMedian = rocCE->GetMedian(rocCEOut);
1516 Double_t pulserMedian = rocPulser->GetMedian(rocCEOut);
1517 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1518 if (rocPulserOut->GetValue(ichannel)>0) {
1519 rocPulser->SetValue(ichannel,pulserMedian);
1520 rocOut->SetValue(ichannel,1);
1522 if (rocCEOut->GetValue(ichannel)>0) {
1523 rocCE->SetValue(ichannel,ceMedian);
1524 rocOut->SetValue(ichannel,1);
1529 // remove pulser time 0
1531 padCE->Add(padPulser,-1);
1536 Float_t chi2Af,chi2Cf;
1537 padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1541 AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1542 padCEFitGY->SetName("padCEFitGy");
1544 AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1545 padCEFit->SetName("padCEFit");
1547 AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE);
1548 padCEDiff->SetName("padCEDiff");
1549 padCEDiff->Add(padCEFit,-1.);
1552 padCE->Add(padCEFitGY,-1.);
1554 padCE->Add(padPulser,1.);
1555 Double_t padmedian = padCE->GetMedian();
1556 padCE->Add(-padmedian); // normalize to median
1558 // Replace outliers by fit value - median of diff per given chamber -GY fit
1560 for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1561 AliTPCCalROC * rocOut = padOut->GetCalROC(iroc);
1562 AliTPCCalROC * rocCE = padCE->GetCalROC(iroc);
1563 AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc);
1564 AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc);
1565 AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc);
1567 Double_t diffMedian = rocCEDiff->GetMedian(rocOut);
1568 for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1569 if (rocOut->GetValue(ichannel)==0) continue;
1570 Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1571 rocCE->SetValue(ichannel,value);
1577 //dump to the file - result can be visualized
1578 AliTPCPreprocessorOnline preprocesor;
1579 preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1580 preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1581 preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1582 preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1584 preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1585 preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1587 preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1588 preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1589 preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1590 preprocesor.DumpToFile(dumpfile);
1593 delete padPulserOut;
1606 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1608 // find the closest point to xref in x direction
1609 // return dx and value
1611 index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1612 if (index<0) index=0;
1613 if (index>=graph->GetN()-1) index=graph->GetN()-2;
1614 if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1615 dx = xref-graph->GetX()[index];
1616 y = graph->GetY()[index];
1621 Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1623 // Get the correction of the trigger offset
1624 // combining information from the laser track calibration
1625 // and from cosmic calibration
1628 // timeStamp - tim stamp in seconds
1629 // deltaT - integration period to calculate offset
1630 // deltaTLaser -max validity of laser data
1631 // valType - 0 - median, 1- mean
1633 // Integration vaues are just recomendation - if not possible to get points
1634 // automatically increase the validity by factor 2
1635 // (recursive algorithm until one month of data taking)
1638 const Float_t kLaserCut=0.0005;
1639 const Int_t kMaxPeriod=3600*24*30*12; // one year max
1640 const Int_t kMinPoints=20;
1642 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1644 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1646 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1647 if (!array) return 0;
1649 TGraphErrors *laserA[3]={0,0,0};
1650 TGraphErrors *laserC[3]={0,0,0};
1651 TGraphErrors *cosmicAll=0;
1652 laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1653 laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1654 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1657 if (!cosmicAll) return 0;
1658 Int_t nmeasC=cosmicAll->GetN();
1659 Float_t *tdelta = new Float_t[nmeasC];
1661 for (Int_t i=0;i<nmeasC;i++){
1662 if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1663 Float_t ccosmic=cosmicAll->GetY()[i];
1664 Double_t yA=0,yC=0,dA=0,dC=0;
1665 if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1666 if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1667 //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1668 //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1670 if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1672 if (TMath::Abs(yA-yC)<kLaserCut) {
1675 if (i%2==0) claser=yA;
1676 if (i%2==1) claser=yC;
1678 tdelta[nused]=ccosmic-claser;
1681 if (nused<kMinPoints &&deltaT<kMaxPeriod) return AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1682 if (nused<kMinPoints) {
1683 printf("AliFatal: No time offset calibration available\n");
1686 Double_t median = TMath::Median(nused,tdelta);
1687 Double_t mean = TMath::Mean(nused,tdelta);
1689 return (valType==0) ? median:mean;
1692 Double_t AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1694 // Get the correction of the drift velocity
1695 // combining information from the laser track calibration
1696 // and from cosmic calibration
1698 // dist - return value - distance to closest point in graph
1700 // timeStamp - tim stamp in seconds
1701 // deltaT - integration period to calculate time0 offset
1702 // deltaTLaser -max validity of laser data
1703 // valType - 0 - median, 1- mean
1705 // Integration vaues are just recomendation - if not possible to get points
1706 // automatically increase the validity by factor 2
1707 // (recursive algorithm until one month of data taking)
1711 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1713 AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1715 array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1716 if (!array) return 0;
1717 TGraphErrors *cosmicAll=0;
1718 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1719 if (!cosmicAll) return 0;
1721 AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1723 Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1724 Double_t vcosmic= AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1725 if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1726 if (timeStamp<cosmicAll->GetX()[0]) vcosmic=cosmicAll->GetY()[0];
1733 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1734 cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1735 laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1737 Double_t *yvd= new Double_t[cosmicAll->GetN()];
1738 Double_t *yt0= new Double_t[cosmicAll->GetN()];
1739 for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1740 for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1742 TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1743 TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1749 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1752 // Create a default name for the gui file
1755 return Form("guiRefTreeRun%s.root",GetRefValidity());
1758 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1761 // Create a gui reference tree
1762 // if dirname and filename are empty default values will be used
1763 // this is the recommended way of using this function
1764 // it allows to check whether a file with the given run validity alredy exists
1766 if (!AliCDBManager::Instance()->GetDefaultStorage()){
1767 AliError("Default Storage not set. Cannot create reference calibration Tree!");
1771 TString file=filename;
1772 if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1774 AliTPCPreprocessorOnline prep;
1775 //noise and pedestals
1776 if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1777 if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1778 if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1780 if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1781 if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1782 if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1783 if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1785 if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1786 if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1787 if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1788 if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1790 if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1791 if (fRefALTROZsThr ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr )));
1792 if (fRefALTROFPED ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED )));
1793 if (fRefALTROAcqStop ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop )));
1794 if (fRefALTROMasked ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked )));
1796 AliTPCdataQA *dataQA=fRefDataQA;
1798 if (dataQA->GetNLocalMaxima())
1799 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1800 if (dataQA->GetMaxCharge())
1801 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1802 if (dataQA->GetMeanCharge())
1803 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1804 if (dataQA->GetNoThreshold())
1805 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1806 if (dataQA->GetNTimeBins())
1807 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1808 if (dataQA->GetNPads())
1809 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1810 if (dataQA->GetTimePosition())
1811 prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1813 prep.DumpToFile(file.Data());
1817 Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1819 // Get the correction of the drift velocity using the laser tracks calbration
1822 // timeStamp - tim stamp in seconds
1823 // deltaT - integration period to calculate time0 offset
1824 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1825 // Note in case no data form both A and C side - the value from active side used
1826 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1827 TGraphErrors *grlaserA=0;
1828 TGraphErrors *grlaserC=0;
1829 Double_t vlaserA=0, vlaserC=0;
1830 if (!array) return 0;
1831 grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1832 grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1835 AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1836 if (TMath::Abs(dist)>deltaT) vlaserA= deltaY;
1837 else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1840 AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1841 if (TMath::Abs(dist)>deltaT) vlaserC= deltaY;
1842 else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1844 if (side==0) return vlaserA;
1845 if (side==1) return vlaserC;
1846 Double_t mdrift=(vlaserA+vlaserC)*0.5;
1847 if (!grlaserA) return vlaserC;
1848 if (!grlaserC) return vlaserA;
1853 Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1855 // Get the correction of the drift velocity using the CE laser data
1856 // combining information from the CE, laser track calibration
1857 // and P/T calibration
1860 // timeStamp - tim stamp in seconds
1861 // deltaT - integration period to calculate time0 offset
1862 // side - 0 - A side, 1 - C side, 2 - mean from both sides
1863 TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime();
1864 if (!arrT) return 0;
1865 AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters();
1866 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
1867 AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1870 Double_t corrPTA = 0, corrPTC=0;
1871 Double_t ltime0A = 0, ltime0C=0;
1873 Double_t corrA=0, corrC=0;
1874 Double_t timeA=0, timeC=0;
1875 TGraph *graphA = (TGraph*)arrT->At(72);
1876 TGraph *graphC = (TGraph*)arrT->At(73);
1877 if (!graphA && !graphC) return 0.;
1878 if (graphA &&graphA->GetN()>0) {
1879 AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1880 timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1881 Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1882 ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1883 if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0);
1884 corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1887 if (graphC&&graphC->GetN()>0){
1888 AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1889 timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1890 Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1891 ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1892 if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0);
1893 corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1897 if (side ==0 ) return corrA;
1898 if (side ==1 ) return corrC;
1899 Double_t corrM= (corrA+corrC)*0.5;
1900 if (!graphA) corrM=corrC;
1901 if (!graphC) corrM=corrA;
1905 Double_t AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1907 // return drift velocity using the TPC-ITS matchin method
1908 // return also distance to the closest point
1910 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1911 TGraphErrors *graph=0;
1913 if (!array) return 0;
1914 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
1915 if (!graph) return 0;
1917 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1918 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1922 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1924 // Get time dependent time 0 (trigger delay in cm) correction
1926 // timestamp - timestamp
1929 // Notice - Extrapolation outside of calibration range - using constant function
1931 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1932 TGraphErrors *graph=0;
1934 if (!array) return 0;
1935 graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
1936 if (!graph) return 0;
1938 AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1939 Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1947 Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1949 // VERY obscure method - we need something in framework
1950 // Find the TPC runs with temperature OCDB entry
1951 // cache the start and end of the run
1953 AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1954 if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1955 if (!storage) return 0;
1956 TString path=storage->GetURI();
1960 if (path.Contains("local")){ // find the list if local system
1961 path.ReplaceAll("local://","");
1962 path+="TPC/Calib/Temperature";
1963 command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data());
1965 runsT=gSystem->GetFromPipe(command);
1967 TObjArray *arr= runsT.Tokenize("r");
1970 TArrayI indexes(arr->GetEntries());
1971 TArrayI runs(arr->GetEntries());
1973 {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1974 Int_t irunN = atoi(arr->At(irun)->GetName());
1975 if (irunN<startRun) continue;
1976 if (irunN>stopRun) continue;
1977 runs[naccept]=irunN;
1981 fRunsStart.Set(fRuns.fN);
1982 fRunsStop.Set(fRuns.fN);
1983 TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1984 for (Int_t irun=0; irun<fRuns.fN; irun++) fRuns[irun]=runs[indexes[irun]];
1987 AliCDBEntry * entry = 0;
1988 {for (Int_t irun=0;irun<fRuns.fN; irun++){
1989 entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1990 if (!entry) continue;
1991 AliTPCSensorTempArray * tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1992 if (!tmpRun) continue;
1993 fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1994 fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1995 // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2001 Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){
2003 // binary search - find the run for given time stamp
2005 Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2006 Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2008 for (Int_t index=index0; index<=index1; index++){
2009 if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2011 printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2014 if (cindex<0) cindex =(index0+index1)/2;
2018 return fRuns[cindex];
2025 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2027 // filter outlyer measurement
2028 // Only points around median +- sigmaCut filtered
2029 if (!graph) return 0;
2031 Int_t npoints0 = graph->GetN();
2034 Double_t *outx=new Double_t[npoints0];
2035 Double_t *outy=new Double_t[npoints0];
2038 if (npoints0<kMinPoints) return 0;
2039 for (Int_t iter=0; iter<3; iter++){
2041 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2042 if (graph->GetY()[ipoint]==0) continue;
2043 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2044 outx[npoints] = graph->GetX()[ipoint];
2045 outy[npoints] = graph->GetY()[ipoint];
2048 if (npoints<=1) break;
2049 medianY =TMath::Median(npoints,outy);
2050 rmsY =TMath::RMS(npoints,outy);
2053 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2058 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2060 // filter outlyer measurement
2061 // Only points around median +- cut filtered
2062 if (!graph) return 0;
2064 Int_t npoints0 = graph->GetN();
2067 Double_t *outx=new Double_t[npoints0];
2068 Double_t *outy=new Double_t[npoints0];
2071 if (npoints0<kMinPoints) return 0;
2072 for (Int_t iter=0; iter<3; iter++){
2074 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2075 if (graph->GetY()[ipoint]==0) continue;
2076 if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2077 outx[npoints] = graph->GetX()[ipoint];
2078 outy[npoints] = graph->GetY()[ipoint];
2081 if (npoints<=1) break;
2082 medianY =TMath::Median(npoints,outy);
2083 rmsY =TMath::RMS(npoints,outy);
2086 if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2092 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2094 // filter outlyer measurement
2095 // Only points with normalized errors median +- sigmaCut filtered
2097 Int_t kMinPoints=10;
2098 Int_t npoints0 = graph->GetN();
2100 Float_t medianErr=0, rmsErr=0;
2101 Double_t *outx=new Double_t[npoints0];
2102 Double_t *outy=new Double_t[npoints0];
2103 Double_t *erry=new Double_t[npoints0];
2104 Double_t *nerry=new Double_t[npoints0];
2105 Double_t *errx=new Double_t[npoints0];
2108 if (npoints0<kMinPoints) return 0;
2109 for (Int_t iter=0; iter<3; iter++){
2111 for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2112 nerry[npoints] = graph->GetErrorY(ipoint);
2113 if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2114 erry[npoints] = graph->GetErrorY(ipoint);
2115 outx[npoints] = graph->GetX()[ipoint];
2116 outy[npoints] = graph->GetY()[ipoint];
2117 errx[npoints] = graph->GetErrorY(ipoint);
2120 if (npoints==0) break;
2121 medianErr=TMath::Median(npoints,erry);
2122 medianY =TMath::Median(npoints,outy);
2123 rmsErr =TMath::RMS(npoints,erry);
2125 TGraphErrors *graphOut=0;
2126 if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2135 void AliTPCcalibDButil::Sort(TGraph *graph){
2137 // sort array - neccessay for approx
2139 Int_t npoints = graph->GetN();
2140 Int_t *indexes=new Int_t[npoints];
2141 Double_t *outx=new Double_t[npoints];
2142 Double_t *outy=new Double_t[npoints];
2143 TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2144 for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2145 for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2146 for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2147 for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2150 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2152 // smmoth graph - mean on the interval
2155 Int_t npoints = graph->GetN();
2156 Double_t *outy=new Double_t[npoints];
2158 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2159 Double_t lx=graph->GetX()[ipoint];
2160 Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2161 Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2162 if (index0<0) index0=0;
2163 if (index1>=npoints-1) index1=npoints-1;
2164 if ((index1-index0)>1){
2165 outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2167 outy[ipoint]=graph->GetY()[ipoint];
2170 // TLinearFitter fitter(3,"pol2");
2171 // for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2172 // Double_t lx=graph->GetX()[ipoint];
2173 // Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2174 // Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2175 // if (index0<0) index0=0;
2176 // if (index1>=npoints-1) index1=npoints-1;
2177 // fitter.ClearPoints();
2178 // for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2179 // if ((index1-index0)>1){
2180 // outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2182 // outy[ipoint]=graph->GetY()[ipoint];
2188 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2189 graph->GetY()[ipoint] = outy[ipoint];
2194 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2196 // Use constant interpolation outside of range
2199 printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2202 if (graph->GetN()<1){
2203 printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2206 if (xref<graph->GetX()[0]) return graph->GetY()[0];
2207 if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2208 return graph->Eval( xref);
2211 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy, Double_t sigmaCut){
2213 // Filter DCS sensor information
2214 // ymin - minimal value
2216 // maxdy - maximal deirivative
2217 // sigmaCut - cut on values and derivative in terms of RMS distribution
2218 // Return value - accepted fraction
2222 // 0. Calculate median and rms of values in specified range
2223 // 1. Filter out outliers - median+-sigmaCut*rms
2224 // values replaced by median
2226 AliSplineFit * fit = sensor->GetFit();
2227 if (!fit) return 0.;
2228 Int_t nknots = fit->GetKnots();
2235 Double_t *yin0 = new Double_t[nknots];
2236 Double_t *yin1 = new Double_t[nknots];
2239 for (Int_t iknot=0; iknot< nknots; iknot++){
2240 if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2241 yin0[naccept] = fit->GetY0()[iknot];
2242 yin1[naccept] = fit->GetY1()[iknot];
2243 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2253 Double_t medianY0=0, medianY1=0;
2254 Double_t rmsY0 =0, rmsY1=0;
2255 medianY0 = TMath::Median(naccept, yin0);
2256 medianY1 = TMath::Median(naccept, yin1);
2257 rmsY0 = TMath::RMS(naccept, yin0);
2258 rmsY1 = TMath::RMS(naccept, yin1);
2261 // 1. Filter out outliers - median+-sigmaCut*rms
2262 // values replaced by median
2263 // if replaced the derivative set to 0
2265 for (Int_t iknot=0; iknot< nknots; iknot++){
2267 if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2268 if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2269 if (nknots<2) fit->GetY1()[iknot]=0;
2270 if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2272 fit->GetY0()[iknot]=medianY0;
2273 fit->GetY1()[iknot]=0;
2280 return Float_t(naccept)/Float_t(nknots);
2283 Float_t AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2285 // Filter temperature array
2286 // tempArray - array of temperatures -
2287 // ymin - minimal accepted temperature - default 15
2288 // ymax - maximal accepted temperature - default 22
2289 // sigmaCut - values filtered on interval median+-sigmaCut*rms - defaut 5
2290 // return value - fraction of filtered sensors
2291 const Double_t kMaxDy=0.1;
2292 Int_t nsensors=tempArray->NumSensors();
2293 if (nsensors==0) return 0.;
2295 for (Int_t isensor=0; isensor<nsensors; isensor++){
2296 AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2297 if (!sensor) continue;
2298 //printf("%d\n",isensor);
2299 FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2300 if (sensor->GetFit()==0){
2302 tempArray->RemoveSensorNum(isensor);
2307 return Float_t(naccept)/Float_t(nsensors);
2311 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2314 // Input parameters:
2315 // deltaT - smoothing window (in seconds)
2316 // cutAbs - max distance of the time info to the median (in time bins)
2317 // cutSigma - max distance (in the RMS)
2318 // pcstream - optional debug streamer to store original and filtered info
2319 // Hardwired parameters:
2320 // kMinPoints =10; // minimal number of points to define the CE
2321 // kMinSectors=12; // minimal number of sectors to define sideCE
2323 // 0. Filter almost emty graphs (kMinPoints=10)
2324 // 1. calculate median and RMS per side
2325 // 2. Filter graphs - in respect with side medians
2326 // - cutAbs and cutDelta used
2327 // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2328 // 4. Calculate mean for A side and C side
2330 const Int_t kMinPoints =10; // minimal number of points to define the CE
2331 const Int_t kMinSectors=12; // minimal number of sectors to define sideCE
2332 const Int_t kMinTime =400; // minimal arrival time of CE
2333 TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2335 TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData();
2336 if (!cearray) return;
2341 AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2342 AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2343 if ( tempMapCE && cavernPressureCE){
2345 // Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2346 // FilterSensor(cavernPressureCE,960,1050,10, 5.);
2347 // if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2350 // recalculate P/T correction map for time of the CE
2351 AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2352 driftCalib->SetName("driftPTCE");
2353 driftCalib->SetTitle("driftPTCE");
2354 cearray->AddLast(driftCalib);
2358 // 0. Filter almost emty graphs
2361 for (Int_t i=0; i<72;i++){
2362 TGraph *graph= (TGraph*)arrT->At(i);
2363 if (!graph) continue;
2364 if (graph->GetN()<kMinPoints){
2366 delete graph; // delete empty graph
2369 if (tmin<0) tmin = graph->GetX()[0];
2370 if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2372 if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2373 if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2376 // 1. calculate median and RMS per side
2378 TArrayF arrA(100000), arrC(100000);
2380 Double_t medianA=0, medianC=0;
2381 Double_t rmsA=0, rmsC=0;
2382 for (Int_t isec=0; isec<72;isec++){
2383 TGraph *graph= (TGraph*)arrT->At(isec);
2384 if (!graph) continue;
2385 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2386 if (graph->GetY()[ipoint]<kMinTime) continue;
2387 if (nA>=arrA.fN) arrA.Set(nA*2);
2388 if (nC>=arrC.fN) arrC.Set(nC*2);
2389 if (isec%36<18) arrA[nA++]= graph->GetY()[ipoint];
2390 if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2394 medianA=TMath::Median(nA,arrA.fArray);
2395 rmsA =TMath::RMS(nA,arrA.fArray);
2398 medianC=TMath::Median(nC,arrC.fArray);
2399 rmsC =TMath::RMS(nC,arrC.fArray);
2402 // 2. Filter graphs - in respect with side medians
2404 TArrayD vecX(100000), vecY(100000);
2405 for (Int_t isec=0; isec<72;isec++){
2406 TGraph *graph= (TGraph*)arrT->At(isec);
2407 if (!graph) continue;
2408 Double_t median = (isec%36<18) ? medianA: medianC;
2409 Double_t rms = (isec%36<18) ? rmsA: rmsC;
2411 for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2412 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2413 if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2414 vecX[naccept]= graph->GetX()[ipoint];
2415 vecY[naccept]= graph->GetY()[ipoint];
2418 if (naccept<kMinPoints){
2419 arrT->AddAt(0,isec);
2420 delete graph; // delete empty graph
2423 TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2425 arrT->AddAt(graph2,isec);
2428 // 3. Cut in respect wit the graph median
2430 for (Int_t i=0; i<72;i++){
2431 TGraph *graph= (TGraph*)arrT->At(i);
2432 if (!graph) continue;
2436 TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2437 if (!graphTS0) continue;
2438 if (graphTS0->GetN()<kMinPoints) {
2444 TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2446 AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2448 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2449 (*pcstream)<<"filterCE"<<
2454 "graphTS0.="<<graphTS0<<
2455 "graphTS.="<<graphTS<<
2459 if (!graphTS) continue;
2460 arrT->AddAt(graphTS,i);
2464 // Recalculate the mean time A side C side
2466 TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2467 Int_t meanPoints=(nA+nC)/72; // mean number of points
2468 for (Int_t itime=0; itime<200; itime++){
2470 Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2471 for (Int_t i=0; i<72;i++){
2472 TGraph *graph= (TGraph*)arrT->At(i);
2473 if (!graph) continue;
2474 if (graph->GetN()<(meanPoints/4)) continue;
2475 if ( (i%36)<18 ) arrA[nA++]=graph->Eval(time);
2476 if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2480 yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2481 yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2482 eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2483 eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2486 Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2487 Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2489 Int_t run = AliTPCcalibDB::Instance()->GetRun();
2490 (*pcstream)<<"filterAC"<<
2499 TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2500 TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2501 TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2502 if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2503 TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2504 if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2507 if (nA<kMinSectors) arrT->AddAt(0,72);
2508 else arrT->AddAt(graphTSA,72);
2509 if (nC<kMinSectors) arrT->AddAt(0,73);
2510 else arrT->AddAt(graphTSC,73);
2514 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2516 // Filter Drift velocity measurement using the tracks
2517 // 0. remove outlyers - error based
2521 const Int_t kMinPoints=1; // minimal number of points to define value
2522 TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2525 for (Int_t i=0; i<arrT->GetEntries();i++){
2526 TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2527 if (!graph) continue;
2528 if (graph->GetN()<kMinPoints){
2533 TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2535 delete graph; arrT->AddAt(0,i); continue;
2537 if (graph2->GetN()<1) {
2538 delete graph; arrT->AddAt(0,i); continue;
2540 graph2->SetName(graph->GetName());
2541 graph2->SetTitle(graph->GetTitle());
2542 arrT->AddAt(graph2,i);
2544 (*pcstream)<<"filterTracks"<<
2549 "graph2.="<<graph2<<
2560 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2563 // get laser time offset
2564 // median around timeStamp+-deltaT
2565 // QA - chi2 needed for later usage - to be added
2566 // - currently cut on error
2569 Double_t kMinDelay=0.01;
2570 Double_t kMinDelayErr=0.0001;
2572 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2573 if (!array) return 0;
2574 TGraphErrors *tlaser=0;
2576 if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2577 if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2579 if (!tlaser) return 0;
2580 Int_t npoints0= tlaser->GetN();
2581 if (npoints0==0) return 0;
2582 Double_t *xlaser = new Double_t[npoints0];
2583 Double_t *ylaser = new Double_t[npoints0];
2585 for (Int_t i=0;i<npoints0;i++){
2587 if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2588 if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2589 xlaser[npoints]=tlaser->GetX()[npoints];
2590 ylaser[npoints]=tlaser->GetY()[npoints];
2595 Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2596 Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2597 //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2598 if (index0<0) index0=0;
2599 if (index1>=npoints-1) index1=npoints-1;
2600 if (index1-index0<kMinPoints) return 0;
2602 //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2603 Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2612 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2614 // Filter Goofie data
2615 // goofieArray - points will be filtered
2616 // deltaT - smmothing time window
2617 // cutSigma - outler sigma cut in rms
2618 // minVn, maxVd- range absolute cut for variable vd/pt
2621 // Ignore goofie if not enough points
2623 const Int_t kMinPoints = 3;
2626 TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2627 TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2628 TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2629 TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2630 if (!graphvd) return;
2631 if (graphvd->GetN()<kMinPoints){
2633 goofieArray->GetSensorNum(2)->SetGraph(0);
2637 // 1. Caluclate medians of critical variables
2643 Double_t medianpt=0;
2644 Double_t medianvd=0, sigmavd=0;
2645 Double_t medianan=0;
2646 Double_t medianaf=0;
2647 Int_t entries=graphvd->GetN();
2648 Double_t yvdn[10000];
2651 for (Int_t ipoint=0; ipoint<entries; ipoint++){
2652 if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2653 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2654 if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2655 yvdn[nvd++]=graphvd->GetY()[ipoint];
2657 if (nvd<kMinPoints){
2659 goofieArray->GetSensorNum(2)->SetGraph(0);
2663 Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2664 if (nuni>=kMinPoints){
2665 AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2667 medianvd = TMath::Median(nvd, yvdn);
2670 TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2671 TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2672 TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2673 TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2674 TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2675 TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2683 // 2. Make outlyer graph
2686 TGraph graphOut(*graphvd);
2687 for (Int_t i=0; i<entries;i++){
2689 Bool_t isOut=kFALSE;
2690 if (graphpt->GetY()[i]<=0.0000001) { graphOut.GetY()[i]=1; continue;}
2691 if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) { graphOut.GetY()[i]=1; continue;}
2693 if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2695 if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2696 if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2697 if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2698 graphOut.GetY()[i]= (isOut)?1:0;
2701 if (nOK<kMinPoints) {
2703 goofieArray->GetSensorNum(2)->SetGraph(0);
2707 // 3. Filter out outlyers - and smooth
2709 TVectorF vmedianArray(goofieArray->NumSensors());
2710 TVectorF vrmsArray(goofieArray->NumSensors());
2711 Double_t xnew[10000];
2712 Double_t ynew[10000];
2714 junk.SetOwner(kTRUE);
2718 for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2720 AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2721 TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2723 if (!sensor) continue;
2724 graphOld = sensor->GetGraph();
2726 sensor->SetGraph(0);
2728 for (Int_t i=0;i<entries;i++){
2729 if (graphOut.GetY()[i]>0.5) continue;
2730 xnew[nused]=graphOld->GetX()[i];
2731 ynew[nused]=graphOld->GetY()[i];
2734 graphNew = new TGraph(nused,xnew,ynew);
2735 junk.AddLast(graphNew);
2736 junk.AddLast(graphOld);
2738 graphNew0 = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2740 junk.AddLast(graphNew0);
2741 graphNew1 = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2743 junk.AddLast(graphNew1);
2744 graphNew2 = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2746 vrmsArray[isensor] =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2747 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2748 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2749 AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2750 printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2751 vmedianArray[isensor]=median;
2757 if (!graphOld) { isOK=kFALSE; graphOld =&graphOut;}
2758 if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2759 if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2760 if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2761 (*pcstream)<<"goofieA"<<
2762 Form("isOK_%d.=",isensor)<<isOK<<
2763 Form("s_%d.=",isensor)<<sensor<<
2764 Form("gr_%d.=",isensor)<<graphOld<<
2765 Form("gr0_%d.=",isensor)<<graphNew0<<
2766 Form("gr1_%d.=",isensor)<<graphNew1<<
2767 Form("gr2_%d.=",isensor)<<graphNew2;
2768 if (isOK) sensor->SetGraph(graphNew2);
2770 (*pcstream)<<"goofieA"<<
2771 "vmed.="<<&vmedianArray<<
2772 "vrms.="<<&vrmsArray<<
2774 junk.Delete(); // delete temoprary graphs
2782 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2784 // Make a statistic matrix
2785 // Input parameters:
2786 // array - TObjArray of AliRelKalmanAlign
2787 // minFraction - minimal ration of accepted tracks
2788 // minStat - minimal statistic (number of accepted tracks)
2789 // maxvd - maximal deviation for the 1
2791 // columns - Mean, Median, RMS
2792 // row - parameter type (rotation[3], translation[3], drift[3])
2793 if (!array) return 0;
2794 if (array->GetEntries()<=0) return 0;
2795 // Int_t entries = array->GetEntries();
2796 Int_t entriesFast = array->GetEntriesFast();
2798 TVectorD *valArray[9];
2799 for (Int_t i=0; i<9; i++){
2800 valArray[i] = new TVectorD(entriesFast);
2803 for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
2804 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
2805 if (!kalman) continue;
2806 if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
2807 if (kalman->GetNUpdates()<minStat) continue;
2808 if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
2809 kalman->GetState(state);
2810 for (Int_t ipar=0; ipar<9; ipar++)
2811 (*valArray[ipar])[naccept]=state[ipar];
2814 if (naccept<2) return 0;
2815 TMatrixD *pstat=new TMatrixD(9,3);
2816 TMatrixD &stat=*pstat;
2817 for (Int_t ipar=0; ipar<9; ipar++){
2818 stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
2819 stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
2820 stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
2826 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
2828 // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
2830 // array - input array
2831 // stat - mean parameters statistic
2833 // sigmaCut - maximal allowed deviation from mean in terms of RMS
2834 if (!array) return 0;
2835 if (array->GetEntries()<=0) return 0;
2836 if (!(&stat)) return 0;
2837 // error increase in 1 hour
2838 const Double_t kerrsTime[9]={
2839 0.00001, 0.00001, 0.00001,
2840 0.001, 0.001, 0.001,
2841 0.0001, 0.001, 0.0001};
2844 Int_t entries = array->GetEntriesFast();
2845 TObjArray *sArray= new TObjArray(entries);
2846 AliRelAlignerKalman * sKalman =0;
2848 for (Int_t i=0; i<entries; i++){
2849 Int_t index=(direction)? entries-i-1:i;
2850 AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
2851 if (!kalman) continue;
2853 kalman->GetState(state);
2854 for (Int_t ipar=0; ipar<9; ipar++){
2855 if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
2857 if (!sKalman &&isOK) {
2858 sKalman=new AliRelAlignerKalman(*kalman);
2859 sKalman->SetRejectOutliers(kFALSE);
2860 sKalman->SetRunNumber(kalman->GetRunNumber());
2861 sKalman->SetTimeStamp(kalman->GetTimeStamp());
2863 if (!sKalman) continue;
2864 Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
2865 for (Int_t ipar=0; ipar<9; ipar++){
2866 // (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
2867 // (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
2868 // (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
2869 (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
2871 sKalman->SetRunNumber(kalman->GetRunNumber());
2872 if (!isOK) sKalman->SetRunNumber(0);
2873 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2874 if (!isOK) continue;
2875 sKalman->SetRejectOutliers(kFALSE);
2876 sKalman->SetRunNumber(kalman->GetRunNumber());
2877 sKalman->SetTimeStamp(kalman->GetTimeStamp());
2878 sKalman->Merge(kalman);
2879 sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2885 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
2887 // Merge 2 RelKalman arrays
2889 // arrayP - rel kalman in direction plus
2890 // arrayM - rel kalman in direction minus
2891 if (!arrayP) return 0;
2892 if (arrayP->GetEntries()<=0) return 0;
2893 if (!arrayM) return 0;
2894 if (arrayM->GetEntries()<=0) return 0;
2895 Int_t entries = arrayP->GetEntriesFast();
2896 TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
2897 for (Int_t i=0; i<entries; i++){
2898 AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
2899 AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
2900 if (!kalmanP) continue;
2901 if (!kalmanM) continue;
2902 AliRelAlignerKalman *kalman = new AliRelAlignerKalman(*kalmanP);
2903 kalman->Merge(kalmanM);
2904 array->AddAt(kalman,i);