/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ /////////////////////////////////////////////////////////////////////////////// // // // Class providing the calculation of derived quantities (mean,rms,fits,...) // // of calibration entries // /* */ //////////////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include #include #include #include "AliTPCcalibDB.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" #include "AliTPCROC.h" #include "AliTPCmapper.h" #include "AliTPCParam.h" #include "AliTPCCalibRaw.h" #include "TGraphErrors.h" #include "AliLog.h" #include "AliTPCcalibDButil.h" #include "AliTPCPreprocessorOnline.h" #include "AliTPCCalibVdrift.h" ClassImp(AliTPCcalibDButil) AliTPCcalibDButil::AliTPCcalibDButil() : TObject(), fCalibDB(AliTPCcalibDB::Instance()), fPadNoise(0x0), fPedestals(0x0), fPulserTmean(0x0), fPulserTrms(0x0), fPulserQmean(0x0), fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")), fCETmean(0x0), fCETrms(0x0), fCEQmean(0x0), fALTROMasked(0x0), fCalibRaw(0x0), fRefPadNoise(0x0), fRefPedestals(0x0), fRefPulserTmean(0x0), fRefPulserTrms(0x0), fRefPulserQmean(0x0), fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")), fRefCETmean(0x0), fRefCETrms(0x0), fRefCEQmean(0x0), fRefALTROMasked(0x0), fRefCalibRaw(0x0), fGoofieArray(0x0), fMapper(new AliTPCmapper(0x0)), fNpulserOutliers(-1), fIrocTimeOffset(0), fCETmaxLimitAbs(1.5), fPulTmaxLimitAbs(1.5), fPulQmaxLimitAbs(5), fPulQminLimit(11), fRuns(0), // run list with OCDB info fRunsStart(0), // start time for given run fRunsStop(0) // stop time for given run { // // Default ctor // } //_____________________________________________________________________________________ AliTPCcalibDButil::~AliTPCcalibDButil() { // // dtor // delete fPulserOutlier; delete fRefPulserOutlier; delete fMapper; if (fRefPadNoise) delete fRefPadNoise; if (fRefPedestals) delete fRefPedestals; if (fRefPulserTmean) delete fRefPulserTmean; if (fRefPulserTrms) delete fRefPulserTrms; if (fRefPulserQmean) delete fRefPulserQmean; if (fRefCETmean) delete fRefCETmean; if (fRefCETrms) delete fRefCETrms; if (fRefCEQmean) delete fRefCEQmean; if (fRefALTROMasked) delete fRefALTROMasked; if (fRefCalibRaw) delete fRefCalibRaw; } //_____________________________________________________________________________________ void AliTPCcalibDButil::UpdateFromCalibDB() { // // Update pointers from calibDB // fPadNoise=fCalibDB->GetPadNoise(); fPedestals=fCalibDB->GetPedestals(); fPulserTmean=fCalibDB->GetPulserTmean(); fPulserTrms=fCalibDB->GetPulserTrms(); fPulserQmean=fCalibDB->GetPulserQmean(); fCETmean=fCalibDB->GetCETmean(); fCETrms=fCalibDB->GetCETrms(); fCEQmean=fCalibDB->GetCEQmean(); fALTROMasked=fCalibDB->GetALTROMasked(); fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun()); fCalibRaw=fCalibDB->GetCalibRaw(); UpdatePulserOutlierMap(); } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC, Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE) { // // Process the CE data for this run // the return TVectorD arrays contian the results of the fit // noutliersCE contains the number of pads marked as outliers, // not including masked and edge pads // //retrieve CE and ALTRO data if (!fCETmean){ TString fitString(fitFormula); fitString.ReplaceAll("++","#"); Int_t ndim=fitString.CountChar('#')+2; fitResultsA.ResizeTo(ndim); fitResultsC.ResizeTo(ndim); fitResultsA.Zero(); fitResultsC.Zero(); noutliersCE=-1; return; } noutliersCE=0; //create outlier map AliTPCCalPad *out=0; if (outCE) out=outCE; else out=new AliTPCCalPad("outCE","outCE"); AliTPCCalROC *rocMasked=0x0; //loop over all channels for (UInt_t iroc=0;irockNsec;++iroc){ AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc); if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc); AliTPCCalROC *rocOut=out->GetCalROC(iroc); if (!rocData) { noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc); rocOut->Add(1.); continue; } //add time offset to IROCs if (irocGetNInnerSector()) rocData->Add(fIrocTimeOffset); //select outliers UInt_t nrows=rocData->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadSetValue(irow,ipad,0); //exclude masked pads if (rocMasked && rocMasked->GetValue(irow,ipad)) { rocOut->SetValue(irow,ipad,1); continue; } //exclude first two rows in IROC and last two rows in OROC if (iroc<36){ if (irow<2) rocOut->SetValue(irow,ipad,1); } else { if (irow>nrows-3) rocOut->SetValue(irow,ipad,1); } //exclude edge pads if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1); Float_t valTmean=rocData->GetValue(irow,ipad); //exclude values that are exactly 0 if (valTmean==0) { rocOut->SetValue(irow,ipad,1); ++noutliersCE; } // exclude channels with too large variations if (TMath::Abs(valTmean)>fCETmaxLimitAbs) { rocOut->SetValue(irow,ipad,1); ++noutliersCE; } } } } //perform fit TMatrixD dummy; Float_t chi2Af,chi2Cf; fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf); chi2A=chi2Af; chi2C=chi2Cf; if (!outCE) delete out; } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian, TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian, Float_t &driftTimeA, Float_t &driftTimeC ) { // // Calculate statistical information from the CE graphs for drift time and charge // //reset arrays vecTEntries.ResizeTo(72); vecTMean.ResizeTo(72); vecTRMS.ResizeTo(72); vecTMedian.ResizeTo(72); vecQEntries.ResizeTo(72); vecQMean.ResizeTo(72); vecQRMS.ResizeTo(72); vecQMedian.ResizeTo(72); vecTEntries.Zero(); vecTMean.Zero(); vecTRMS.Zero(); vecTMedian.Zero(); vecQEntries.Zero(); vecQMean.Zero(); vecQRMS.Zero(); vecQMedian.Zero(); driftTimeA=0; driftTimeC=0; TObjArray *arrT=fCalibDB->GetCErocTtime(); TObjArray *arrQ=fCalibDB->GetCErocQtime(); if (arrT){ for (Int_t isec=0;isec<74;++isec){ TGraph *gr=(TGraph*)arrT->At(isec); if (!gr) continue; TVectorD values; Int_t npoints = gr->GetN(); values.ResizeTo(npoints); Int_t nused =0; //skip first points, theres always some problems with finding the CE position for (Int_t ipoint=4; ipointGetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){ values[nused]=gr->GetY()[ipoint]; nused++; } } // if (isec<72) vecTEntries[isec]= nused; if (nused>1){ if (isec<72){ vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray()); vecTMean[isec] = TMath::Mean(nused,values.GetMatrixArray()); vecTRMS[isec] = TMath::RMS(nused,values.GetMatrixArray()); } else if (isec==72){ driftTimeA=TMath::Median(nused,values.GetMatrixArray()); } else if (isec==73){ driftTimeC=TMath::Median(nused,values.GetMatrixArray()); } } } } if (arrQ){ for (Int_t isec=0;isecGetEntriesFast();++isec){ TGraph *gr=(TGraph*)arrQ->At(isec); if (!gr) continue; TVectorD values; Int_t npoints = gr->GetN(); values.ResizeTo(npoints); Int_t nused =0; for (Int_t ipoint=0; ipointGetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){ values[nused]=gr->GetY()[ipoint]; nused++; } } // vecQEntries[isec]= nused; if (nused>1){ vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray()); vecQMean[isec] = TMath::Mean(nused,values.GetMatrixArray()); vecQRMS[isec] = TMath::RMS(nused,values.GetMatrixArray()); } } } } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions, TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions, Int_t &nonMaskedZero) { // // process noise data // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1], // OROCs small pads [2] and OROCs large pads [3] // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot) // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error // //set proper size and reset const UInt_t infoSize=4; vNoiseMean.ResizeTo(infoSize); vNoiseMeanSenRegions.ResizeTo(infoSize); vNoiseRMS.ResizeTo(infoSize); vNoiseRMSSenRegions.ResizeTo(infoSize); vNoiseMean.Zero(); vNoiseMeanSenRegions.Zero(); vNoiseRMS.Zero(); vNoiseRMSSenRegions.Zero(); nonMaskedZero=0; //counters TVectorD c(infoSize); TVectorD cs(infoSize); //tpc parameters AliTPCParam par; par.Update(); //retrieve noise and ALTRO data if (!fPadNoise) return; AliTPCCalROC *rocMasked=0x0; //create IROC, OROC1, OROC2 and sensitive region masks for (UInt_t isec=0;isecGetCalROC(isec); if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec); UInt_t nrows=noiseROC->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad)) continue; Float_t noiseVal=noiseROC->GetValue(irow,ipad); //check if noise==0 if (noiseVal==0) { ++nonMaskedZero; continue; } //check for nan if ( !(noiseVal<10000000) ){ printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal); continue; } Int_t cpad=(Int_t)ipad-(Int_t)npads/2; Int_t masksen=1; // sensitive pards are not masked (0) if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive) if (isecGetNInnerSector()){ //IROCs if (irow>19&&irow<46){ if (TMath::Abs(cpad)<7) masksen=0; //IROC spot } Int_t type=1; vNoiseMean[type]+=noiseVal; vNoiseRMS[type]+=noiseVal*noiseVal; ++c[type]; if (!masksen){ vNoiseMeanSenRegions[type]+=noiseVal; vNoiseRMSSenRegions[type]+=noiseVal*noiseVal; ++cs[type]; } } else { //OROCs //define sensive regions if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive if ( irow>75 ){ Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad); if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive } if ((Int_t)irowverySmall){ // printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]); mean=vNoiseMean[i]/c[i]; rms=vNoiseRMS[i]; rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean)); } vNoiseMean[i]=mean; vNoiseRMS[i]=rms; if (cs[i]>verySmall){ meanSen=vNoiseMeanSenRegions[i]/cs[i]; rmsSen=vNoiseRMSSenRegions[i]; rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen)); } vNoiseMeanSenRegions[i]=meanSen; vNoiseRMSSenRegions[i]=rmsSen; } } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime) { // // Process the Pulser information // vMeanTime: pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C // const UInt_t infoSize=4; //reset counters to error number vMeanTime.ResizeTo(infoSize); vMeanTime.Zero(); //counter TVectorD c(infoSize); //retrieve pulser and ALTRO data if (!fPulserTmean) return; // //get Outliers AliTPCCalROC *rocOut=0x0; for (UInt_t isec=0;isecGetCalROC(isec); if (!tmeanROC) continue; rocOut=fPulserOutlier->GetCalROC(isec); UInt_t nchannels=tmeanROC->GetNchannels(); for (UInt_t ichannel=0;ichannelGetValue(ichannel)) continue; Float_t val=tmeanROC->GetValue(ichannel); Int_t type=isec/18; vMeanTime[type]+=val; ++c[type]; } } //calculate mean for (UInt_t itype=0; itype0) vMeanTime[itype]/=c[itype]; else vMeanTime[itype]=0; } } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked) { // // Get Values from ALTRO configuration data // nMasked=-1; if (!fALTROMasked) return; nMasked=0; for (Int_t isec=0;iseckNsec; ++isec){ AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec); for (UInt_t ichannel=0; ichannelGetNchannels();++ichannel){ if (rocMasked->GetValue(ichannel)) ++nMasked; } } } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS) { // // Proces Goofie values, return statistical information of the currently set goofieArray // The meaning of the entries are given below /* 1 TPC_ANODE_I_A00_STAT 2 TPC_DVM_CO2 3 TPC_DVM_DriftVelocity 4 TPC_DVM_FCageHV 5 TPC_DVM_GainFar 6 TPC_DVM_GainNear 7 TPC_DVM_N2 8 TPC_DVM_NumberOfSparks 9 TPC_DVM_PeakAreaFar 10 TPC_DVM_PeakAreaNear 11 TPC_DVM_PeakPosFar 12 TPC_DVM_PeakPosNear 13 TPC_DVM_PickupHV 14 TPC_DVM_Pressure 15 TPC_DVM_T1_Over_P 16 TPC_DVM_T2_Over_P 17 TPC_DVM_T_Over_P 18 TPC_DVM_TemperatureS1 */ if (!fGoofieArray){ Int_t nsensors=19; vecEntries.ResizeTo(nsensors); vecMedian.ResizeTo(nsensors); vecMean.ResizeTo(nsensors); vecRMS.ResizeTo(nsensors); vecEntries.Zero(); vecMedian.Zero(); vecMean.Zero(); vecRMS.Zero(); return; } Double_t kEpsilon=0.0000000001; Double_t kBig=100000000000.; Int_t nsensors = fGoofieArray->NumSensors(); vecEntries.ResizeTo(nsensors); vecMedian.ResizeTo(nsensors); vecMean.ResizeTo(nsensors); vecRMS.ResizeTo(nsensors); TVectorF values; for (Int_t isensor=0; isensorNumSensors();isensor++){ AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor); if (gsensor && gsensor->GetGraph()){ Int_t npoints = gsensor->GetGraph()->GetN(); // filter zeroes values.ResizeTo(npoints); Int_t nused =0; for (Int_t ipoint=0; ipointGetGraph()->GetY()[ipoint])>kEpsilon && TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])GetGraph()->GetY()[ipoint]; nused++; } } // vecEntries[isensor]= nused; if (nused>1){ vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray()); vecMean[isensor] = TMath::Mean(nused,values.GetMatrixArray()); vecRMS[isensor] = TMath::RMS(nused,values.GetMatrixArray()); } } } } //_____________________________________________________________________________________ void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations) { // // check the variations of the pedestal data to the reference pedestal data // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively. // const Int_t npar=4; TVectorF vThres(npar); //thresholds Int_t nActive=0; //number of active channels //reset and set thresholds pedestalDeviations.ResizeTo(npar); for (Int_t i=0;iGetCalROC(isec); AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec); AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec); AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec); UInt_t nrows=mROC->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad)) continue; if (mRefROC->GetValue(irow,ipad)) continue; Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad)); for (Int_t i=0;ivThres[i]) ++pedestalDeviations.GetMatrixArray()[i]; } ++nActive; }//end ipad }//ind irow }//end isec if (nActive>0){ for (Int_t i=0;iGetCalROC(isec); AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec); AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec); AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec); UInt_t nrows=mROC->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad)) continue; if (mRefROC->GetValue(irow,ipad)) continue; Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1; for (Int_t i=0;ivThres[i]) ++noiseDeviations.GetMatrixArray()[i]; } ++nActive; }//end ipad }//ind irow }//end isec if (nActive>0){ for (Int_t i=0;iGetCalROC(isec); AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec); AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec); // AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec); AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec); AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec); AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec); Float_t pt_mean=ptROC->GetMean(oROC); UInt_t nrows=mROC->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad)) continue; if (mRefROC->GetValue(irow,ipad)) continue; //don't user edge pads if (ipad==0||ipad==npads-1) continue; //data Float_t pq=pqROC->GetValue(irow,ipad); Float_t pqRef=pqRefROC->GetValue(irow,ipad); Float_t pt=ptROC->GetValue(irow,ipad); // Float_t ptRef=ptRefROC->GetValue(irow,ipad); //comparisons q Float_t deviation=TMath::Abs(pq/pqRef-1); for (Int_t i=0;ivThres[i]) ++pulserQdeviations.GetMatrixArray()[i]; } if (pqRef>11&&pq<11) ++npadsOffAdd; varQMean+=pq-pqRef; //comparisons t if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB; ++nActive; }//end ipad }//ind irow }//end isec if (nActive>0){ for (Int_t i=0;iMultiply(0.); fNpulserOutliers=-1; return; } AliTPCCalROC *rocMasked=0x0; fNpulserOutliers=0; //Create Outlier Map for (UInt_t isec=0;isecGetCalROC(isec); AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec); AliTPCCalROC *outROC=pulOut->GetCalROC(isec); if (!tmeanROC||!qmeanROC) { //reset outliers in this ROC outROC->Multiply(0.); continue; } if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec); // Double_t dummy=0; // Float_t qmedian=qmeanROC->GetLTM(&dummy,.5); // Float_t tmedian=tmeanROC->GetLTM(&dummy,.5); UInt_t nrows=tmeanROC->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad); Float_t t=tmeanROC->GetValue(irow,ipad); //masked channels are outliers if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1; //edge pads are outliers if (ipad==0||ipad==npads-1) masked=1; //channels with too large charge or timing deviation from the meadian are outliers // if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1; if (qSetValue(irow,ipad,outlier+masked); fNpulserOutliers+=outlier; } } } } //_____________________________________________________________________________________ AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C ) { // // Create pad time0 object from pulser and/or CE data, depending on the selected model // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser // Model 2: use CE data and a combination CE fit + pulser in the outlier regions. // // In case model 2 is invoked - gy arival time gradient is also returned // gyA=0; gyC=0; AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model)); // decide between different models if (model==0||model==1){ TVectorD vMean; if (model==1) ProcessPulser(vMean); for (UInt_t isec=0;isecGetCalROC(isec); if (!rocPulTmean) continue; AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec); AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec); Float_t mean=rocPulTmean->GetMean(rocOut); //treat case where a whole partition is masked if (mean==0) mean=rocPulTmean->GetMean(); if (model==1) { Int_t type=isec/18; mean=vMean[type]; } UInt_t nrows=rocTime0->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad); //in case of an outlier pad use the mean of the altro values. //This should be the most precise guess in that case. if (rocOut->GetValue(irow,ipad)) { time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut); if (time==0) time=mean; } Float_t val=time-mean; rocTime0->SetValue(irow,ipad,val); } } } } else if (model==2){ Double_t pgya,pgyc,pchi2a,pchi2c; AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c); fCETmean->Add(padPulser,-1.); TVectorD vA,vC; AliTPCCalPad outCE("outCE","outCE"); Int_t nOut; ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE); AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC); // AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC); if (!padFit) { delete padPulser; return 0;} gyA=vA[2]; gyC=vC[2]; fCETmean->Add(padPulser,1.); padTime0->Add(fCETmean); padTime0->Add(padFit,-1); delete padPulser; TVectorD vFitROC; TMatrixD mFitROC; Float_t chi2; for (UInt_t isec=0;isecGetCalROC(isec); AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec); AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec); AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec); rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2); AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec); Float_t mean=rocPulTmean->GetMean(rocOutPul); if (mean==0) mean=rocPulTmean->GetMean(); UInt_t nrows=rocTime0->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadGetValue(irow,ipad)-mean; if (rocOutCE->GetValue(irow,ipad)){ Float_t valOut=rocCEfit->GetValue(irow,ipad); if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser; rocTime0->SetValue(irow,ipad,valOut); } } } delete rocCEfit; } delete padFit; } Double_t median = padTime0->GetMedian(); padTime0->Add(-median); // normalize to median return padTime0; } //_____________________________________________________________________________________ Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut) { if (roc==0) return 0.; const Int_t sector=roc->GetSector(); AliTPCROC *tpcRoc=AliTPCROC::Instance(); const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad); Float_t mean=0; Int_t n=0; //loop over a small range around the requested pad (+-10 rows/pads) for (Int_t irow=row-10;irow(Int_t)tpcRoc->GetNRows(sector)-1) continue; for (Int_t ipad=pad-10; ipad(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue; const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad); if (altroRoc!=altroCurr) continue; if ( rocOut && rocOut->GetValue(irow,ipad) ) continue; Float_t val=roc->GetValue(irow,ipad); mean+=val; ++n; } } if (n>0) mean/=n; return mean; } //_____________________________________________________________________________________ void AliTPCcalibDButil::SetRefFile(const char* filename) { // // load cal pad objects form the reference file // TDirectory *currDir=gDirectory; TFile f(filename); fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals"); fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise"); //pulser data fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean"); fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms"); fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean"); //CE data fRefCETmean=(AliTPCCalPad*)f.Get("CETmean"); fRefCETrms=(AliTPCCalPad*)f.Get("CETrms"); fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean"); //Altro data // fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart"); // fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr"); // fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED"); // fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop"); fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked"); f.Close(); currDir->cd(); } AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin, Float_t cutTrmsMax, Float_t cutMaxDistT){ // // Author: marian.ivanov@cern.ch // // Create outlier map for CE study // Parameters: // Return value - outlyer map // noutlyersCE - number of outlyers // minSignal - minimal total Q signal // cutRMSMin - minimal width of the signal in respect to the median // cutRMSMax - maximal width of the signal in respect to the median // cutMaxDistT - maximal deviation from time median per chamber // // Outlyers criteria: // 0. Exclude masked pads // 1. Exclude first two rows in IROC and last two rows in OROC // 2. Exclude edge pads // 3. Exclude channels with too large variations // 4. Exclude pads with too small signal // 5. Exclude signal with outlyers RMS // 6. Exclude channels to far from the chamber median noutliersCE=0; //create outlier map AliTPCCalPad *out=ceOut; if (!out) out= new AliTPCCalPad("outCE","outCE"); AliTPCCalROC *rocMasked=0x0; if (!fCETmean) return 0; if (!fCETrms) return 0; if (!fCEQmean) return 0; // //loop over all channels // Double_t rmsMedian = fCETrms->GetMedian(); for (UInt_t iroc=0;irockNsec;++iroc){ AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc); if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc); AliTPCCalROC *rocOut = out->GetCalROC(iroc); AliTPCCalROC *rocCEQ = fCEQmean->GetCalROC(iroc); AliTPCCalROC *rocCETrms = fCETrms->GetCalROC(iroc); Double_t trocMedian = rocData->GetMedian(); // if (!rocData) { noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc); rocOut->Add(1.); continue; } // //select outliers UInt_t nrows=rocData->GetNrows(); for (UInt_t irow=0;irowGetNPads(irow); for (UInt_t ipad=0;ipadSetValue(irow,ipad,0); Float_t valTmean=rocData->GetValue(irow,ipad); Float_t valQmean=rocCEQ->GetValue(irow,ipad); Float_t valTrms =rocCETrms->GetValue(irow,ipad); //0. exclude masked pads if (rocMasked && rocMasked->GetValue(irow,ipad)) { rocOut->SetValue(irow,ipad,1); continue; } //1. exclude first two rows in IROC and last two rows in OROC if (iroc<36){ if (irow<2) rocOut->SetValue(irow,ipad,1); } else { if (irow>nrows-3) rocOut->SetValue(irow,ipad,1); } //2. exclude edge pads if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1); //exclude values that are exactly 0 if (valTmean==0) { rocOut->SetValue(irow,ipad,1); ++noutliersCE; } //3. exclude channels with too large variations if (TMath::Abs(valTmean)>fCETmaxLimitAbs) { rocOut->SetValue(irow,ipad,1); ++noutliersCE; } // //4. exclude channels with too small signal if (valQmeanSetValue(irow,ipad,1); ++noutliersCE; } // //5. exclude channels with too small rms if (valTrmscutTrmsMax*rmsMedian){ rocOut->SetValue(irow,ipad,1); ++noutliersCE; } // //6. exclude channels to far from the chamber median if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){ rocOut->SetValue(irow,ipad,1); ++noutliersCE; } } } } // return out; } AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){ // // Author: marian.ivanov@cern.ch // // Create outlier map for Pulser // Parameters: // Return value - outlyer map // noutlyersPulser - number of outlyers // cutTime - absolute cut - distance to the median of chamber // cutnRMSQ - nsigma cut from median q distribution per chamber // cutnRMSrms - nsigma cut from median rms distribution // Outlyers criteria: // 0. Exclude masked pads // 1. Exclude time outlyers (default 3 time bins) // 2. Exclude q outlyers (default 5 sigma) // 3. Exclude rms outlyers (default 5 sigma) noutliersPulser=0; AliTPCCalPad *out=pulserOut; if (!out) out= new AliTPCCalPad("outPulser","outPulser"); AliTPCCalROC *rocMasked=0x0; if (!fPulserTmean) return 0; if (!fPulserTrms) return 0; if (!fPulserQmean) return 0; // //loop over all channels // for (UInt_t iroc=0;irockNsec;++iroc){ if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc); AliTPCCalROC *rocData = fPulserTmean->GetCalROC(iroc); AliTPCCalROC *rocOut = out->GetCalROC(iroc); AliTPCCalROC *rocPulserQ = fPulserQmean->GetCalROC(iroc); AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc); // Double_t rocMedianT = rocData->GetMedian(); Double_t rocMedianQ = rocPulserQ->GetMedian(); Double_t rocRMSQ = rocPulserQ->GetRMS(); Double_t rocMedianTrms = rocPulserTrms->GetMedian(); Double_t rocRMSTrms = rocPulserTrms->GetRMS(); for (UInt_t ichannel=0;ichannelGetNchannels();++ichannel){ rocOut->SetValue(ichannel,0); Float_t valTmean=rocData->GetValue(ichannel); Float_t valQmean=rocPulserQ->GetValue(ichannel); Float_t valTrms =rocPulserTrms->GetValue(ichannel); Int_t isOut=0; if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1; if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1; if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1; rocOut->SetValue(ichannel,isOut); if (isOut) noutliersPulser++; } } return out; } AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){ // // Author : Marian Ivanov // Create pad time0 correction map using information from the CE and from pulser // // // Return PadTime0 to be used for time0 relative alignment // if dump file specified intermediat results are dumped to the fiel and can be visualized // using $ALICE_ROOT/TPC/script/gui application // // fitResultsA - fitParameters A side // fitResultsC - fitParameters C side // chi2A - chi2/ndf for A side (assuming error 1 time bin) // chi2C - chi2/ndf for C side (assuming error 1 time bin) // // // Algorithm: // 1. Find outlier map for CE // 2. Find outlier map for Pulser // 3. Replace outlier by median at given sector (median without outliers) // 4. Substract from the CE data pulser // 5. Fit the CE with formula // 5.1) (IROC-OROC) offset // 5.2) gx // 5.3) gy // 5.4) (lx-xmid) // 5.5) (IROC-OROC)*(lx-xmid) // 5.6) (ly/lx)^2 // 6. Substract gy fit dependence from the CE data // 7. Add pulser back to CE data // 8. Replace outliers by fit value - median of diff per given chamber -GY fit // 9. return CE data // // Time0 <= padCE = padCEin -padCEfitGy - if not outlier // Time0 <= padCE = padFitAll-padCEfitGy - if outlier // fit formula 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)"; // output for fit formula 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)"; // gy part of formula 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)"; // // if (!fCETmean) return 0; Double_t pgya,pgyc,pchi2a,pchi2c; AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut); AliTPCCalPad * padCEOut = CreateCEOutlyerMap(nOut); AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c); AliTPCCalPad * padCE = new AliTPCCalPad(*fCETmean); AliTPCCalPad * padCEIn = new AliTPCCalPad(*fCETmean); AliTPCCalPad * padOut = new AliTPCCalPad("padOut","padOut"); padPulser->SetName("padPulser"); padPulserOut->SetName("padPulserOut"); padCE->SetName("padCE"); padCEIn->SetName("padCEIn"); padCEOut->SetName("padCEOut"); padOut->SetName("padOut"); // // make combined outlyers map // and replace outlyers in maps with median for chamber // for (UInt_t iroc=0;irockNsec;++iroc){ AliTPCCalROC * rocOut = padOut->GetCalROC(iroc); AliTPCCalROC * rocPulser = padPulser->GetCalROC(iroc); AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc); AliTPCCalROC * rocCEOut = padCEOut->GetCalROC(iroc); AliTPCCalROC * rocCE = padCE->GetCalROC(iroc); Double_t ceMedian = rocCE->GetMedian(rocCEOut); Double_t pulserMedian = rocPulser->GetMedian(rocCEOut); for (UInt_t ichannel=0;ichannelGetNchannels();++ichannel){ if (rocPulserOut->GetValue(ichannel)>0) { rocPulser->SetValue(ichannel,pulserMedian); rocOut->SetValue(ichannel,1); } if (rocCEOut->GetValue(ichannel)>0) { rocCE->SetValue(ichannel,ceMedian); rocOut->SetValue(ichannel,1); } } } // // remove pulser time 0 // padCE->Add(padPulser,-1); // // Make fits // TMatrixD dummy; Float_t chi2Af,chi2Cf; padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf); chi2A=chi2Af; chi2C=chi2Cf; // AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC); padCEFitGY->SetName("padCEFitGy"); // AliTPCCalPad *padCEFit =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC); padCEFit->SetName("padCEFit"); // AliTPCCalPad* padCEDiff = new AliTPCCalPad(*padCE); padCEDiff->SetName("padCEDiff"); padCEDiff->Add(padCEFit,-1.); // // padCE->Add(padCEFitGY,-1.); padCE->Add(padPulser,1.); Double_t padmedian = padCE->GetMedian(); padCE->Add(-padmedian); // normalize to median // // Replace outliers by fit value - median of diff per given chamber -GY fit // for (UInt_t iroc=0;irockNsec;++iroc){ AliTPCCalROC * rocOut = padOut->GetCalROC(iroc); AliTPCCalROC * rocCE = padCE->GetCalROC(iroc); AliTPCCalROC * rocCEFit = padCEFit->GetCalROC(iroc); AliTPCCalROC * rocCEFitGY = padCEFitGY->GetCalROC(iroc); AliTPCCalROC * rocCEDiff = padCEDiff->GetCalROC(iroc); // Double_t diffMedian = rocCEDiff->GetMedian(rocOut); for (UInt_t ichannel=0;ichannelGetNchannels();++ichannel){ if (rocOut->GetValue(ichannel)==0) continue; Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian; rocCE->SetValue(ichannel,value); } } // // if (dumpfile){ //dump to the file - result can be visualized AliTPCPreprocessorOnline preprocesor; preprocesor.AddComponent(new AliTPCCalPad(*padCE)); preprocesor.AddComponent(new AliTPCCalPad(*padCEIn)); preprocesor.AddComponent(new AliTPCCalPad(*padCEFit)); preprocesor.AddComponent(new AliTPCCalPad(*padOut)); // preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY)); preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff)); // preprocesor.AddComponent(new AliTPCCalPad(*padCEOut)); preprocesor.AddComponent(new AliTPCCalPad(*padPulser)); preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut)); preprocesor.DumpToFile(dumpfile); } delete padPulser; delete padPulserOut; delete padCEIn; delete padCEOut; delete padOut; delete padCEDiff; delete padCEFitGY; return padCE; } Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){ // // find the closest point to xref in x direction // return dx and value Int_t index=0; index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref); if (index<0) index=0; if (index>=graph->GetN()-1) index=graph->GetN()-2; if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++; dx = xref-graph->GetX()[index]; y = graph->GetY()[index]; return index; } Double_t AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){ // // Get the correction of the trigger offset // combining information from the laser track calibration // and from cosmic calibration // // run - run number // timeStamp - tim stamp in seconds // deltaT - integration period to calculate offset // deltaTLaser -max validity of laser data // valType - 0 - median, 1- mean // // Integration vaues are just recomendation - if not possible to get points // automatically increase the validity by factor 2 // (recursive algorithm until one month of data taking) // // const Float_t kLaserCut=0.0005; const Int_t kMaxPeriod=3600*24*30*3; // 3 month max const Int_t kMinPoints=20; // TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); if (!array) { AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE); } array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); if (!array) return 0; // TGraphErrors *laserA[3]={0,0,0}; TGraphErrors *laserC[3]={0,0,0}; TGraphErrors *cosmicAll=0; laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A"); laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C"); cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL"); // // if (!cosmicAll) return 0; Int_t nmeasC=cosmicAll->GetN(); Float_t *tdelta = new Float_t[nmeasC]; Int_t nused=0; for (Int_t i=0;iGetX()[i]-timeStamp)>deltaT) continue; Float_t ccosmic=cosmicAll->GetY()[i]; Double_t yA=0,yC=0,dA=0,dC=0; if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA); if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC); //yA=laserA[1]->Eval(cosmicAll->GetX()[i]); //yC=laserC[1]->Eval(cosmicAll->GetX()[i]); // if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue; Float_t claser=0; if (TMath::Abs(yA-yC)GetTimeVdriftSplineRun(run); if (!array) { AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE); } array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); if (!array) return 0; TGraphErrors *cosmicAll=0; cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL"); if (!cosmicAll) return 0; Double_t grY=0; AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY); Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType); Double_t vcosmic= AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp); if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1]) vcosmic=timeStamp>cosmicAll->GetY()[cosmicAll->GetN()-1]; if (timeStampGetX()[0]) vcosmic=cosmicAll->GetY()[0]; return vcosmic+t0; /* Example usage: Int_t run=89000 TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL"); laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A"); // Double_t *yvd= new Double_t[cosmicAll->GetN()]; Double_t *yt0= new Double_t[cosmicAll->GetN()]; for (Int_t i=0; iGetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]); for (Int_t i=0; iGetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]); TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd); TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0); */ } Double_t AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){ // // Get the correction of the drift velocity using the laser tracks calbration // // run - run number // timeStamp - tim stamp in seconds // deltaT - integration period to calculate time0 offset // side - 0 - A side, 1 - C side, 2 - mean from both sides // Note in case no data form both A and C side - the value from active side used TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); TGraphErrors *grlaserA=0; TGraphErrors *grlaserC=0; Double_t vlaserA=0, vlaserC=0; if (!array) return 0; grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A"); grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C"); Double_t deltaY; if (grlaserA) { AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY); if (TMath::Abs(dist)>deltaT) vlaserA= deltaY; else vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp); } if (grlaserC) { AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY); if (TMath::Abs(dist)>deltaT) vlaserC= deltaY; else vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp); } if (side==0) return vlaserA; if (side==1) return vlaserC; Double_t mdrift=(vlaserA+vlaserC)*0.5; if (!grlaserA) return vlaserC; if (!grlaserC) return vlaserA; return mdrift; } Double_t AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){ // // Get the correction of the drift velocity using the CE laser data // combining information from the CE, laser track calibration // and P/T calibration // // run - run number // timeStamp - tim stamp in seconds // deltaT - integration period to calculate time0 offset // side - 0 - A side, 1 - C side, 2 - mean from both sides TObjArray *arrT =AliTPCcalibDB::Instance()->GetCErocTtime(); AliTPCParam *param =AliTPCcalibDB::Instance()->GetParameters(); TObjArray* cearray =AliTPCcalibDB::Instance()->GetCEData(); AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE"); // // Double_t corrPTA = 0, corrPTC=0; Double_t ltime0A = 0, ltime0C=0; Double_t gry=0; Double_t corrA=0, corrC=0; Double_t timeA=0, timeC=0; TGraph *graphA = (TGraph*)arrT->At(72); TGraph *graphC = (TGraph*)arrT->At(73); if (!graphA && !graphC) return 0.; if (graphA &&graphA->GetN()>0) { AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry); timeA = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp); Int_t mtime =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5); ltime0A = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0); if (driftCalib) corrPTA = driftCalib->GetPTRelative(timeStamp,0); corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-3.*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1; corrA-=corrPTA; } if (graphC&&graphC->GetN()>0){ AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry); timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp); Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5); ltime0C = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0); if (driftCalib) corrPTC = driftCalib->GetPTRelative(timeStamp,0); corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-3.*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1; corrC-=corrPTC; } if (side ==0 ) return corrA; if (side ==1 ) return corrC; Double_t corrM= (corrA+corrC)*0.5; if (!graphA) corrM=corrC; if (!graphC) corrM=corrA; return corrM; } Int_t AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){ // // VERY obscure method - we need something in framework // Find the TPC runs with temperature OCDB entry // cache the start and end of the run // AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature"); if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage(); if (!storage) return 0; TString path=storage->GetURI(); TString runsT; { TString command; if (path.Contains("local")){ // find the list if local system path.ReplaceAll("local://",""); path+="TPC/Calib/Temperature"; command=Form("ls %s | sed s/_/\\ /g | awk '{print \"r\"$2}' ",path.Data()); } runsT=gSystem->GetFromPipe(command); } TObjArray *arr= runsT.Tokenize("r"); if (!arr) return 0; // TArrayI indexes(arr->GetEntries()); TArrayI runs(arr->GetEntries()); Int_t naccept=0; {for (Int_t irun=0;irunGetEntries();irun++){ Int_t irunN = atoi(arr->At(irun)->GetName()); if (irunNstopRun) continue; runs[naccept]=irunN; naccept++; }} fRuns.Set(naccept); fRunsStart.Set(fRuns.fN); fRunsStop.Set(fRuns.fN); TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE); for (Int_t irun=0; irunGet("TPC/Calib/Temperature",fRuns[irun]); if (!entry) continue; AliTPCSensorTempArray * tmpRun = dynamic_cast(entry->GetObject()); if (!tmpRun) continue; fRunsStart[irun]=tmpRun->GetStartTime().GetSec(); fRunsStop[irun]=tmpRun->GetEndTime().GetSec(); // printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec()); }} return fRuns.fN; } Int_t AliTPCcalibDButil::FindRunTPC(Int_t itime, Bool_t debug){ // // binary search - find the run for given time stamp // Int_t index0 = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime); Int_t index1 = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime); Int_t cindex = -1; for (Int_t index=index0; index<=index1; index++){ if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index; if (debug) { printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime); } } if (cindex<0) cindex =(index0+index1)/2; if (cindex<0) { return 0; } return fRuns[cindex]; } TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){ // // filter outlyer measurement // Only points around median +- sigmaCut filtered if (!graph) return 0; Int_t kMinPoints=2; Int_t npoints0 = graph->GetN(); Int_t npoints=0; Float_t rmsY=0; Double_t *outx=new Double_t[npoints0]; Double_t *outy=new Double_t[npoints0]; // // if (npoints0GetY()[ipoint]==0) continue; if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue; outx[npoints] = graph->GetX()[ipoint]; outy[npoints] = graph->GetY()[ipoint]; npoints++; } if (npoints<=1) break; medianY =TMath::Median(npoints,outy); rmsY =TMath::RMS(npoints,outy); } TGraph *graphOut=0; if (npoints>1) graphOut= new TGraph(npoints,outx,outy); return graphOut; } TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){ // // filter outlyer measurement // Only points around median +- cut filtered if (!graph) return 0; Int_t kMinPoints=2; Int_t npoints0 = graph->GetN(); Int_t npoints=0; Float_t rmsY=0; Double_t *outx=new Double_t[npoints0]; Double_t *outy=new Double_t[npoints0]; // // if (npoints0GetY()[ipoint]==0) continue; if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue; outx[npoints] = graph->GetX()[ipoint]; outy[npoints] = graph->GetY()[ipoint]; npoints++; } if (npoints<=1) break; medianY =TMath::Median(npoints,outy); rmsY =TMath::RMS(npoints,outy); } TGraph *graphOut=0; if (npoints>1) graphOut= new TGraph(npoints,outx,outy); return graphOut; } TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){ // // filter outlyer measurement // Only points with normalized errors median +- sigmaCut filtered // Int_t kMinPoints=10; Int_t npoints0 = graph->GetN(); Int_t npoints=0; Float_t medianErr=0, rmsErr=0; Double_t *outx=new Double_t[npoints0]; Double_t *outy=new Double_t[npoints0]; Double_t *erry=new Double_t[npoints0]; Double_t *nerry=new Double_t[npoints0]; Double_t *errx=new Double_t[npoints0]; // // if (npoints0GetErrorY(ipoint); if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue; erry[npoints] = graph->GetErrorY(ipoint); outx[npoints] = graph->GetX()[ipoint]; outy[npoints] = graph->GetY()[ipoint]; errx[npoints] = graph->GetErrorY(ipoint); npoints++; } if (npoints==0) break; medianErr=TMath::Median(npoints,erry); medianY =TMath::Median(npoints,outy); rmsErr =TMath::RMS(npoints,erry); } TGraphErrors *graphOut=0; if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry); delete []outx; delete []outy; delete []errx; delete []erry; return graphOut; } void AliTPCcalibDButil::Sort(TGraph *graph){ // // sort array - neccessay for approx // Int_t npoints = graph->GetN(); Int_t *indexes=new Int_t[npoints]; Double_t *outx=new Double_t[npoints]; Double_t *outy=new Double_t[npoints]; TMath::Sort(npoints, graph->GetX(),indexes,kFALSE); for (Int_t i=0;iGetX()[indexes[i]]; for (Int_t i=0;iGetY()[indexes[i]]; for (Int_t i=0;iGetX()[i]=outx[i]; for (Int_t i=0;iGetY()[i]=outy[i]; } void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){ // // smmoth graph - mean on the interval // Sort(graph); Int_t npoints = graph->GetN(); Double_t *outy=new Double_t[npoints]; for (Int_t ipoint=0; ipointGetX()[ipoint]; Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta); Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta); if (index0<0) index0=0; if (index1>=npoints-1) index1=npoints-1; if ((index1-index0)>1){ outy[ipoint] = TMath::Mean(index1-index0, &(graph->GetY()[index0])); }else{ outy[ipoint]=graph->GetY()[ipoint]; } } for (Int_t ipoint=0; ipointGetY()[ipoint] = outy[ipoint]; } delete[] outy; } Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){ // // Use constant interpolation outside of range // if (!graph) { printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n"); return 0; } if (graph->GetN()<1){ printf("AliTPCcalibDButil::EvalGraphConst: Empty graph"); return 0; } if (xrefGetX()[0]) return graph->GetY()[0]; if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1]; return graph->Eval( xref); } void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){ // // Filter CE data // 0. remove outlyers // 0.1 absolute cut // 0.2 nsigma cut // 1. smooth the graphs // const Int_t kMinPoints=5; // minimal number of points to define the CE TObjArray *arrT=fCalibDB->GetCErocTtime(); Double_t medianY=0; TObjArray* cearray = fCalibDB->GetCEData(); if (!cearray) return; AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap"); AliDCSSensor * cavernPressureCE = (AliDCSSensor *) cearray->FindObject("CavernPressure"); if ( tempMapCE && cavernPressureCE){ // recalculate P/T correction map for time of the CE AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0); driftCalib->SetName("driftPTCE"); driftCalib->SetTitle("driftPTCE"); cearray->AddLast(driftCalib); } for (Int_t i=0; iGetEntries();i++){ TGraph *graph= (TGraph*)arrT->At(i); if (!graph) continue; if (graph->GetN()AddAt(0,i); delete graph; // delete empty graph continue; } TGraph* graphTS0= AliTPCcalibDButil::FilterGraphMedianAbs(graph,cutAbs,medianY); if (!graphTS0) continue; if (graphTS0->GetN()AddAt(0,i); continue; } TGraph* graphTS= AliTPCcalibDButil::FilterGraphMedian(graphTS0,cutSigma,medianY); graphTS->Sort(); AliTPCcalibDButil::SmoothGraph(graphTS,deltaT); if (pcstream){ Int_t run = AliTPCcalibDB::Instance()->GetRun(); (*pcstream)<<"filterCE"<< "run="<AddAt(graphTS,i); delete graph; } } void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){ // // Filter Drift velocity measurement using the tracks // 0. remove outlyers - error based // cutSigma // // const Int_t kMinPoints=1; // minimal number of points to define value TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); Double_t medianY=0; if (!arrT) return; for (Int_t i=0; iGetEntries();i++){ TGraphErrors *graph= (TGraphErrors*)arrT->At(i); if (!graph) continue; if (graph->GetN()AddAt(0,i); continue; } TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY); if (!graph2) { delete graph; arrT->AddAt(0,i); continue; } if (graph2->GetN()<1) { delete graph; arrT->AddAt(0,i); continue; } graph2->SetName(graph->GetName()); graph2->SetTitle(graph->GetTitle()); arrT->AddAt(graph2,i); if (pcstream){ (*pcstream)<<"filterTracks"<< "run="<NumSensors();isensor++){ // // AliDCSSensor *sensor = goofieArray->GetSensor(isensor); Double_t medianY; if (sensor && sensor->GetGraph()){ TGraph * graph = sensor->GetGraph(); if (isensor==3 && graph->GetN()>1){ // drift velocity TGraph * graphv = AliTPCcalibDButil::FilterGraphMedianAbs(graph,0.2,medianY); delete graph; graph=graphv; } TGraph * graph2 = AliTPCcalibDButil::FilterGraphMedian(graph,cutSigma,medianY); if (!graph2) continue; AliTPCcalibDButil::SmoothGraph(graph2,deltaT); if (pcstream){ Int_t run = AliTPCcalibDB::Instance()->GetRun(); (*pcstream)<<"filterG"<< "run="<GetN()<2) {delete graph2; continue;} delete graph; sensor->SetGraph(graph2); } } } Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){ // // // get laser time offset // median around timeStamp+-deltaT // QA - chi2 needed for later usage - to be added // - currently cut on error // Int_t kMinPoints=1; Double_t kMinDelay=0.01; Double_t kMinDelayErr=0.0001; TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run); if (!array) return 0; TGraphErrors *tlaser=0; if (array){ if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A"); if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C"); } if (!tlaser) return 0; Int_t npoints0= tlaser->GetN(); if (npoints0==0) return 0; Double_t *xlaser = new Double_t[npoints0]; Double_t *ylaser = new Double_t[npoints0]; Int_t npoints=0; for (Int_t i=0;iGetY()[i]<=kMinDelay) continue; // filter zeros if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue; xlaser[npoints]=tlaser->GetX()[npoints]; ylaser[npoints]=tlaser->GetY()[npoints]; npoints++; } // // Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1; Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1; //if (index1-index0 =npoints-1) index1=npoints-1; if (index1-index0