TPC/AliTPCcalibDButil.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   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  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // Class providing the calculation of derived quantities (mean,rms,fits,...) //
  20 //       of calibration entries                                              //
  21 /*
  22
  23
  24 */
  25 ////////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TMath.h>
  28 #include <TVectorT.h>
  29 #include <TObjArray.h>
  30 #include <TGraph.h>
  31 #include <TFile.h>
  32 #include <TDirectory.h>
  33 #include <TMap.h>
  34 #include <TGraphErrors.h>
  35 #include <AliCDBStorage.h>
  36 #include <AliDCSSensorArray.h>
  37 #include <AliTPCSensorTempArray.h>
  38 #include <AliDCSSensor.h>
  39 #include <AliLog.h>
  40 #include <AliCDBEntry.h>
  41 #include <AliCDBManager.h>
  42 #include <AliCDBId.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"
  52 #include "AliLog.h"
  53 #include "AliTPCcalibDButil.h"
  54 #include "AliTPCCalibVdrift.h"
  55 #include "AliMathBase.h"
  56 #include "AliRelAlignerKalman.h"
  57
  58 const Float_t kAlmost0=1.e-30;
  59
  60 ClassImp(AliTPCcalibDButil)
  61 AliTPCcalibDButil::AliTPCcalibDButil() :
  62   TObject(),
  63   fCalibDB(0),
  64   fPadNoise(0x0),
  65   fPedestals(0x0),
  66   fPulserTmean(0x0),
  67   fPulserTrms(0x0),
  68   fPulserQmean(0x0),
  69   fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
  70   fCETmean(0x0),
  71   fCETrms(0x0),
  72   fCEQmean(0x0),
  73   fALTROMasked(0x0),
  74   fCalibRaw(0x0),
  75   fDataQA(0x0),
  76   fRefMap(0x0),
  77   fCurrentRefMap(0x0),
  78   fRefValidity(""),
  79   fRefPadNoise(0x0),
  80   fRefPedestals(0x0),
  81   fRefPedestalMasked(0x0),
  82   fRefPulserTmean(0x0),
  83   fRefPulserTrms(0x0),
  84   fRefPulserQmean(0x0),
  85   fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
  86   fRefPulserMasked(0x0),
  87   fRefCETmean(0x0),
  88   fRefCETrms(0x0),
  89   fRefCEQmean(0x0),
  90   fRefCEMasked(0x0),
  91   fRefALTROFPED(0x0),
  92   fRefALTROZsThr(0x0),
  93   fRefALTROAcqStart(0x0),
  94   fRefALTROAcqStop(0x0),
  95   fRefALTROMasked(0x0),
  96   fRefCalibRaw(0x0),
  97   fRefDataQA(0x0),
  98   fGoofieArray(0x0),
  99   fMapper(new AliTPCmapper(0x0)),
 100   fNpulserOutliers(-1),
 101   fIrocTimeOffset(0),
 102   fCETmaxLimitAbs(1.5),
 103   fPulTmaxLimitAbs(1.5),
 104   fPulQmaxLimitAbs(5),
 105   fPulQminLimit(11),
 106   fRuns(0),                         // run list with OCDB info
 107   fRunsStart(0),                    // start time for given run
 108   fRunsStop(0)                     // stop time for given run
 109 {
 110   //
 111   // Default ctor
 112   //
 113 }
 114 //_____________________________________________________________________________________
 115 AliTPCcalibDButil::~AliTPCcalibDButil()
 116 {
 117   //
 118   // dtor
 119   //
 120   delete fPulserOutlier;
 121   delete fRefPulserOutlier;
 122   delete fMapper;
 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;
 141 }
 142 //_____________________________________________________________________________________
 143 void AliTPCcalibDButil::UpdateFromCalibDB()
 144 {
 145   //
 146   // Update pointers from calibDB
 147   //
 148   if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
 149   fPadNoise=fCalibDB->GetPadNoise();
 150   fPedestals=fCalibDB->GetPedestals();
 151   fPulserTmean=fCalibDB->GetPulserTmean();
 152   fPulserTrms=fCalibDB->GetPulserTrms();
 153   fPulserQmean=fCalibDB->GetPulserQmean();
 154   fCETmean=fCalibDB->GetCETmean();
 155   fCETrms=fCalibDB->GetCETrms();
 156   fCEQmean=fCalibDB->GetCEQmean();
 157   fALTROMasked=fCalibDB->GetALTROMasked();
 158   fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
 159   fCalibRaw=fCalibDB->GetCalibRaw();
 160   fDataQA=fCalibDB->GetDataQA();
 161   UpdatePulserOutlierMap();
 162 //   SetReferenceRun();
 163 //   UpdateRefDataFromOCDB();
 164 }
 165 //_____________________________________________________________________________________
 166 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
 167                                       Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
 168 {
 169   //
 170   // Process the CE data for this run
 171   // the return TVectorD arrays contian the results of the fit
 172   // noutliersCE contains the number of pads marked as outliers,
 173   //   not including masked and edge pads
 174   //
 175
 176   //retrieve CE and ALTRO data
 177   if (!fCETmean){
 178     TString fitString(fitFormula);
 179     fitString.ReplaceAll("++","#");
 180     Int_t ndim=fitString.CountChar('#')+2;
 181     fitResultsA.ResizeTo(ndim);
 182     fitResultsC.ResizeTo(ndim);
 183     fitResultsA.Zero();
 184     fitResultsC.Zero();
 185     noutliersCE=-1;
 186     return;
 187   }
 188   noutliersCE=0;
 189   //create outlier map
 190   AliTPCCalPad *out=0;
 191   if (outCE) out=outCE;
 192   else out=new AliTPCCalPad("outCE","outCE");
 193   AliTPCCalROC *rocMasked=0x0;
 194   //loop over all channels
 195   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 196     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 197     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
 198     AliTPCCalROC *rocOut=out->GetCalROC(iroc);
 199     if (!rocData) {
 200       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
 201       rocOut->Add(1.);
 202       continue;
 203     }
 204     //add time offset to IROCs
 205     if (iroc<AliTPCROC::Instance()->GetNInnerSector())
 206       rocData->Add(fIrocTimeOffset);
 207     //select outliers
 208     UInt_t nrows=rocData->GetNrows();
 209     for (UInt_t irow=0;irow<nrows;++irow){
 210       UInt_t npads=rocData->GetNPads(irow);
 211       for (UInt_t ipad=0;ipad<npads;++ipad){
 212         rocOut->SetValue(irow,ipad,0);
 213         //exclude masked pads
 214         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
 215           rocOut->SetValue(irow,ipad,1);
 216           continue;
 217         }
 218         //exclude first two rows in IROC and last two rows in OROC
 219         if (iroc<36){
 220           if (irow<2) rocOut->SetValue(irow,ipad,1);
 221         } else {
 222           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
 223         }
 224         //exclude edge pads
 225         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
 226         Float_t valTmean=rocData->GetValue(irow,ipad);
 227         //exclude values that are exactly 0
 228         if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
 229           rocOut->SetValue(irow,ipad,1);
 230           ++noutliersCE;
 231         }
 232         // exclude channels with too large variations
 233         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
 234           rocOut->SetValue(irow,ipad,1);
 235           ++noutliersCE;
 236         }
 237       }
 238     }
 239   }
 240   //perform fit
 241   TMatrixD dummy;
 242   Float_t chi2Af,chi2Cf;
 243   fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
 244   chi2A=chi2Af;
 245   chi2C=chi2Cf;
 246   if (!outCE) delete out;
 247 }
 248 //_____________________________________________________________________________________
 249 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
 250                      TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
 251                      Float_t &driftTimeA, Float_t &driftTimeC )
 252 {
 253   //
 254   // Calculate statistical information from the CE graphs for drift time and charge
 255   //
 256
 257   //reset arrays
 258   vecTEntries.ResizeTo(72);
 259   vecTMean.ResizeTo(72);
 260   vecTRMS.ResizeTo(72);
 261   vecTMedian.ResizeTo(72);
 262   vecQEntries.ResizeTo(72);
 263   vecQMean.ResizeTo(72);
 264   vecQRMS.ResizeTo(72);
 265   vecQMedian.ResizeTo(72);
 266   vecTEntries.Zero();
 267   vecTMean.Zero();
 268   vecTRMS.Zero();
 269   vecTMedian.Zero();
 270   vecQEntries.Zero();
 271   vecQMean.Zero();
 272   vecQRMS.Zero();
 273   vecQMedian.Zero();
 274   driftTimeA=0;
 275   driftTimeC=0;
 276   TObjArray *arrT=fCalibDB->GetCErocTtime();
 277   TObjArray *arrQ=fCalibDB->GetCErocQtime();
 278   if (arrT){
 279     for (Int_t isec=0;isec<74;++isec){
 280       TGraph *gr=(TGraph*)arrT->At(isec);
 281       if (!gr) continue;
 282       TVectorD values;
 283       Int_t npoints = gr->GetN();
 284       values.ResizeTo(npoints);
 285       Int_t nused =0;
 286       //skip first points, theres always some problems with finding the CE position
 287       for (Int_t ipoint=4; ipoint<npoints; ipoint++){
 288         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
 289           values[nused]=gr->GetY()[ipoint];
 290           nused++;
 291         }
 292       }
 293       //
 294       if (isec<72) vecTEntries[isec]= nused;
 295       if (nused>1){
 296         if (isec<72){
 297           vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 298           vecTMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 299           vecTRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 300         } else if (isec==72){
 301           driftTimeA=TMath::Median(nused,values.GetMatrixArray());
 302         } else if (isec==73){
 303           driftTimeC=TMath::Median(nused,values.GetMatrixArray());
 304         }
 305       }
 306     }
 307   }
 308   if (arrQ){
 309     for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
 310       TGraph *gr=(TGraph*)arrQ->At(isec);
 311       if (!gr) continue;
 312       TVectorD values;
 313       Int_t npoints = gr->GetN();
 314       values.ResizeTo(npoints);
 315       Int_t nused =0;
 316       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 317         if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
 318           values[nused]=gr->GetY()[ipoint];
 319           nused++;
 320         }
 321       }
 322       //
 323       vecQEntries[isec]= nused;
 324       if (nused>1){
 325         vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 326         vecQMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 327         vecQRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 328       }
 329     }
 330   }
 331 }
 332
 333 //_____________________________________________________________________________________
 334 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
 335                       TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
 336                       Int_t &nonMaskedZero, Int_t &nNaN)
 337 {
 338   //
 339   // process noise data
 340   // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
 341   //    OROCs small pads [2] and OROCs large pads [3]
 342   // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
 343   // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
 344   //
 345
 346   //set proper size and reset
 347   const UInt_t infoSize=4;
 348   vNoiseMean.ResizeTo(infoSize);
 349   vNoiseMeanSenRegions.ResizeTo(infoSize);
 350   vNoiseRMS.ResizeTo(infoSize);
 351   vNoiseRMSSenRegions.ResizeTo(infoSize);
 352   vNoiseMean.Zero();
 353   vNoiseMeanSenRegions.Zero();
 354   vNoiseRMS.Zero();
 355   vNoiseRMSSenRegions.Zero();
 356   nonMaskedZero=0;
 357   nNaN=0;
 358   //counters
 359   TVectorD c(infoSize);
 360   TVectorD cs(infoSize);
 361   //tpc parameters
 362   AliTPCParam par;
 363   par.Update();
 364   //retrieve noise and ALTRO data
 365   if (!fPadNoise) return;
 366   AliTPCCalROC *rocMasked=0x0;
 367   //create IROC, OROC1, OROC2 and sensitive region masks
 368   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 369     AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
 370     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 371     UInt_t nrows=noiseROC->GetNrows();
 372     for (UInt_t irow=0;irow<nrows;++irow){
 373       UInt_t npads=noiseROC->GetNPads(irow);
 374       for (UInt_t ipad=0;ipad<npads;++ipad){
 375         //don't use masked channels;
 376         if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
 377         Float_t noiseVal=noiseROC->GetValue(irow,ipad);
 378         //check if noise==0
 379         if (noiseVal<kAlmost0) {
 380           ++nonMaskedZero;
 381           continue;
 382         }
 383         //check for nan
 384         if ( !(noiseVal<10000000) ){
 385 //           printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
 386           ++nNaN;
 387           continue;
 388         }
 389         Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
 390         Int_t masksen=1; // sensitive pards are not masked (0)
 391         if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
 392         if (isec<AliTPCROC::Instance()->GetNInnerSector()){
 393           //IROCs
 394           if (irow>19&&irow<46){
 395             if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
 396           }
 397           Int_t type=1;
 398           vNoiseMean[type]+=noiseVal;
 399           vNoiseRMS[type]+=noiseVal*noiseVal;
 400           ++c[type];
 401           if (!masksen){
 402             vNoiseMeanSenRegions[type]+=noiseVal;
 403             vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 404             ++cs[type];
 405           }
 406         } else {
 407           //OROCs
 408           //define sensive regions
 409           if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
 410           if ( irow>75 ){
 411             Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
 412             if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
 413           }
 414           if ((Int_t)irow<par.GetNRowUp1()){
 415             //OROC1
 416             Int_t type=2;
 417             vNoiseMean[type]+=noiseVal;
 418             vNoiseRMS[type]+=noiseVal*noiseVal;
 419             ++c[type];
 420             if (!masksen){
 421               vNoiseMeanSenRegions[type]+=noiseVal;
 422               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 423               ++cs[type];
 424             }
 425           }else{
 426             //OROC2
 427             Int_t type=3;
 428             vNoiseMean[type]+=noiseVal;
 429             vNoiseRMS[type]+=noiseVal*noiseVal;
 430             ++c[type];
 431             if (!masksen){
 432               vNoiseMeanSenRegions[type]+=noiseVal;
 433               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 434               ++cs[type];
 435             }
 436           }
 437         }
 438         //whole tpc
 439         Int_t type=0;
 440         vNoiseMean[type]+=noiseVal;
 441         vNoiseRMS[type]+=noiseVal*noiseVal;
 442         ++c[type];
 443         if (!masksen){
 444           vNoiseMeanSenRegions[type]+=noiseVal;
 445           vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 446           ++cs[type];
 447         }
 448       }//end loop pads
 449     }//end loop rows
 450   }//end loop sectors (rocs)
 451
 452   //calculate mean and RMS
 453   const Double_t verySmall=0.0000000001;
 454   for (UInt_t i=0;i<infoSize;++i){
 455     Double_t mean=0;
 456     Double_t rms=0;
 457     Double_t meanSen=0;
 458     Double_t rmsSen=0;
 459
 460     if (c[i]>verySmall){
 461 //       printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
 462       mean=vNoiseMean[i]/c[i];
 463       rms=vNoiseRMS[i];
 464       rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
 465     }
 466     vNoiseMean[i]=mean;
 467     vNoiseRMS[i]=rms;
 468
 469     if (cs[i]>verySmall){
 470       meanSen=vNoiseMeanSenRegions[i]/cs[i];
 471       rmsSen=vNoiseRMSSenRegions[i];
 472       rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
 473     }
 474     vNoiseMeanSenRegions[i]=meanSen;
 475     vNoiseRMSSenRegions[i]=rmsSen;
 476   }
 477 }
 478
 479 //_____________________________________________________________________________________
 480 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
 481 {
 482   //
 483   // Process the Pulser information
 484   // vMeanTime:     pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
 485   //
 486
 487   const UInt_t infoSize=4;
 488   //reset counters to error number
 489   vMeanTime.ResizeTo(infoSize);
 490   vMeanTime.Zero();
 491   //counter
 492   TVectorD c(infoSize);
 493   //retrieve pulser and ALTRO data
 494   if (!fPulserTmean) return;
 495   //
 496   //get Outliers
 497   AliTPCCalROC *rocOut=0x0;
 498   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 499     AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
 500     if (!tmeanROC) continue;
 501     rocOut=fPulserOutlier->GetCalROC(isec);
 502     UInt_t nchannels=tmeanROC->GetNchannels();
 503     for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
 504       if (rocOut && rocOut->GetValue(ichannel)) continue;
 505       Float_t val=tmeanROC->GetValue(ichannel);
 506       Int_t type=isec/18;
 507       vMeanTime[type]+=val;
 508       ++c[type];
 509     }
 510   }
 511   //calculate mean
 512   for (UInt_t itype=0; itype<infoSize; ++itype){
 513     if (c[itype]>0) vMeanTime[itype]/=c[itype];
 514     else vMeanTime[itype]=0;
 515   }
 516 }
 517 //_____________________________________________________________________________________
 518 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
 519 {
 520   //
 521   // Get Values from ALTRO configuration data
 522   //
 523   nMasked=-1;
 524   if (!fALTROMasked) return;
 525   nMasked=0;
 526   for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
 527     AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
 528     for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
 529       if (rocMasked->GetValue(ichannel)) ++nMasked;
 530     }
 531   }
 532 }
 533 //_____________________________________________________________________________________
 534 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
 535 {
 536   //
 537   // Proces Goofie values, return statistical information of the currently set goofieArray
 538   // The meaning of the entries are given below
 539   /*
 540   1       TPC_ANODE_I_A00_STAT
 541   2       TPC_DVM_CO2
 542   3       TPC_DVM_DriftVelocity
 543   4       TPC_DVM_FCageHV
 544   5       TPC_DVM_GainFar
 545   6       TPC_DVM_GainNear
 546   7       TPC_DVM_N2
 547   8       TPC_DVM_NumberOfSparks
 548   9       TPC_DVM_PeakAreaFar
 549   10      TPC_DVM_PeakAreaNear
 550   11      TPC_DVM_PeakPosFar
 551   12      TPC_DVM_PeakPosNear
 552   13      TPC_DVM_PickupHV
 553   14      TPC_DVM_Pressure
 554   15      TPC_DVM_T1_Over_P
 555   16      TPC_DVM_T2_Over_P
 556   17      TPC_DVM_T_Over_P
 557   18      TPC_DVM_TemperatureS1
 558    */
 559   if (!fGoofieArray){
 560     Int_t nsensors=19;
 561     vecEntries.ResizeTo(nsensors);
 562     vecMedian.ResizeTo(nsensors);
 563     vecMean.ResizeTo(nsensors);
 564     vecRMS.ResizeTo(nsensors);
 565     vecEntries.Zero();
 566     vecMedian.Zero();
 567     vecMean.Zero();
 568     vecRMS.Zero();
 569     return;
 570   }
 571   Double_t kEpsilon=0.0000000001;
 572   Double_t kBig=100000000000.;
 573   Int_t nsensors = fGoofieArray->NumSensors();
 574   vecEntries.ResizeTo(nsensors);
 575   vecMedian.ResizeTo(nsensors);
 576   vecMean.ResizeTo(nsensors);
 577   vecRMS.ResizeTo(nsensors);
 578   TVectorF values;
 579   for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
 580     AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
 581     if (gsensor &&  gsensor->GetGraph()){
 582       Int_t npoints = gsensor->GetGraph()->GetN();
 583       // filter zeroes
 584       values.ResizeTo(npoints);
 585       Int_t nused =0;
 586       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 587         if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
 588             TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
 589               values[nused]=gsensor->GetGraph()->GetY()[ipoint];
 590               nused++;
 591             }
 592       }
 593       //
 594       vecEntries[isensor]= nused;
 595       if (nused>1){
 596         vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
 597         vecMean[isensor]   = TMath::Mean(nused,values.GetMatrixArray());
 598         vecRMS[isensor]    = TMath::RMS(nused,values.GetMatrixArray());
 599       }
 600     }
 601   }
 602 }
 603 //_____________________________________________________________________________________
 604 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
 605 {
 606   //
 607   // check the variations of the pedestal data to the reference pedestal data
 608   // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
 609   //
 610   const Int_t npar=4;
 611   TVectorF vThres(npar); //thresholds
 612   Int_t nActive=0;       //number of active channels
 613
 614   //reset and set thresholds
 615   pedestalDeviations.ResizeTo(npar);
 616   for (Int_t i=0;i<npar;++i){
 617     pedestalDeviations.GetMatrixArray()[i]=0;
 618     vThres.GetMatrixArray()[i]=(i+1)*.5;
 619   }
 620   //check all needed data is available
 621   if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
 622   //loop over all channels
 623   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 624     AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
 625     AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
 626     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 627     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 628     UInt_t nrows=mROC->GetNrows();
 629     for (UInt_t irow=0;irow<nrows;++irow){
 630       UInt_t npads=mROC->GetNPads(irow);
 631       for (UInt_t ipad=0;ipad<npads;++ipad){
 632         //don't use masked channels;
 633         if (mROC   ->GetValue(irow,ipad)) continue;
 634         if (mRefROC->GetValue(irow,ipad)) continue;
 635         Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
 636         for (Int_t i=0;i<npar;++i){
 637           if (deviation>vThres[i])
 638             ++pedestalDeviations.GetMatrixArray()[i];
 639         }
 640         ++nActive;
 641       }//end ipad
 642     }//ind irow
 643   }//end isec
 644   if (nActive>0){
 645     for (Int_t i=0;i<npar;++i){
 646       pedestalDeviations.GetMatrixArray()[i]/=nActive;
 647     }
 648   }
 649 }
 650 //_____________________________________________________________________________________
 651 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
 652 {
 653   //
 654   // check the variations of the noise data to the reference noise data
 655   // thresholds are 5, 10, 15 and 20 percent respectively.
 656   //
 657   const Int_t npar=4;
 658   TVectorF vThres(npar); //thresholds
 659   Int_t nActive=0;       //number of active channels
 660
 661   //reset and set thresholds
 662   noiseDeviations.ResizeTo(npar);
 663   for (Int_t i=0;i<npar;++i){
 664     noiseDeviations.GetMatrixArray()[i]=0;
 665     vThres.GetMatrixArray()[i]=(i+1)*.05;
 666   }
 667   //check all needed data is available
 668   if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
 669   //loop over all channels
 670   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 671     AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
 672     AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
 673     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 674     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 675     UInt_t nrows=mROC->GetNrows();
 676     for (UInt_t irow=0;irow<nrows;++irow){
 677       UInt_t npads=mROC->GetNPads(irow);
 678       for (UInt_t ipad=0;ipad<npads;++ipad){
 679         //don't use masked channels;
 680         if (mROC   ->GetValue(irow,ipad)) continue;
 681         if (mRefROC->GetValue(irow,ipad)) continue;
 682         if (nRefROC->GetValue(irow,ipad)==0) continue;
 683         Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
 684         for (Int_t i=0;i<npar;++i){
 685           if (deviation>vThres[i])
 686             ++noiseDeviations.GetMatrixArray()[i];
 687         }
 688         ++nActive;
 689       }//end ipad
 690     }//ind irow
 691   }//end isec
 692   if (nActive>0){
 693     for (Int_t i=0;i<npar;++i){
 694       noiseDeviations.GetMatrixArray()[i]/=nActive;
 695     }
 696   }
 697 }
 698 //_____________________________________________________________________________________
 699 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
 700                                                 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
 701 {
 702   //
 703   // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
 704   // thresholds are .5, 1, 5 and 10 percent respectively.
 705   //
 706   //
 707   const Int_t npar=4;
 708   TVectorF vThres(npar); //thresholds
 709   Int_t nActive=0;       //number of active channels
 710
 711   //reset and set thresholds
 712   pulserQdeviations.ResizeTo(npar);
 713   for (Int_t i=0;i<npar;++i){
 714     pulserQdeviations.GetMatrixArray()[i]=0;
 715   }
 716   npadsOutOneTB=0;
 717   npadsOffAdd=0;
 718   varQMean=0;
 719   vThres.GetMatrixArray()[0]=.005;
 720   vThres.GetMatrixArray()[1]=.01;
 721   vThres.GetMatrixArray()[2]=.05;
 722   vThres.GetMatrixArray()[3]=.1;
 723   //check all needed data is available
 724   if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
 725   //
 726   UpdateRefPulserOutlierMap();
 727   //loop over all channels
 728   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 729     AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
 730     AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
 731     AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
 732 //     AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
 733     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 734     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 735     AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
 736     Float_t ptmean=ptROC->GetMean(oROC);
 737     UInt_t nrows=mROC->GetNrows();
 738     for (UInt_t irow=0;irow<nrows;++irow){
 739       UInt_t npads=mROC->GetNPads(irow);
 740       for (UInt_t ipad=0;ipad<npads;++ipad){
 741         //don't use masked channels;
 742         if (mROC   ->GetValue(irow,ipad)) continue;
 743         if (mRefROC->GetValue(irow,ipad)) continue;
 744         //don't user edge pads
 745         if (ipad==0||ipad==npads-1) continue;
 746         //data
 747         Float_t pq=pqROC->GetValue(irow,ipad);
 748         Float_t pqRef=pqRefROC->GetValue(irow,ipad);
 749         Float_t pt=ptROC->GetValue(irow,ipad);
 750 //         Float_t ptRef=ptRefROC->GetValue(irow,ipad);
 751         //comparisons q
 752         Float_t deviation=TMath::Abs(pq/pqRef-1);
 753         for (Int_t i=0;i<npar;++i){
 754           if (deviation>vThres[i])
 755             ++pulserQdeviations.GetMatrixArray()[i];
 756         }
 757         if (pqRef>11&&pq<11) ++npadsOffAdd;
 758         varQMean+=pq-pqRef;
 759         //comparisons t
 760         if (TMath::Abs(pt-ptmean)>1) ++npadsOutOneTB;
 761         ++nActive;
 762       }//end ipad
 763     }//ind irow
 764   }//end isec
 765   if (nActive>0){
 766     for (Int_t i=0;i<npar;++i){
 767       pulserQdeviations.GetMatrixArray()[i]/=nActive;
 768       varQMean/=nActive;
 769     }
 770   }
 771 }
 772 //_____________________________________________________________________________________
 773 void AliTPCcalibDButil::UpdatePulserOutlierMap()
 774 {
 775   //
 776   // Update the outlier map of the pulser data
 777   //
 778   PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
 779 }
 780 //_____________________________________________________________________________________
 781 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
 782 {
 783   //
 784   // Update the outlier map of the pulser reference data
 785   //
 786   PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
 787 }
 788 //_____________________________________________________________________________________
 789 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
 790 {
 791   //
 792   // Create a map that contains outliers from the Pulser calibration data.
 793   // The outliers include masked channels, edge pads and pads with
 794   //   too large timing and charge variations.
 795   // fNpulserOutliers is the number of outliers in the Pulser calibration data.
 796   //   those do not contain masked and edge pads
 797   //
 798   if (!pulT||!pulQ) {
 799     //reset map
 800     pulOut->Multiply(0.);
 801     fNpulserOutliers=-1;
 802     return;
 803   }
 804   AliTPCCalROC *rocMasked=0x0;
 805   fNpulserOutliers=0;
 806
 807   //Create Outlier Map
 808   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 809     AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
 810     AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
 811     AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
 812     if (!tmeanROC||!qmeanROC) {
 813       //reset outliers in this ROC
 814       outROC->Multiply(0.);
 815       continue;
 816     }
 817     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 818 //     Double_t dummy=0;
 819 //     Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
 820 //     Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
 821     UInt_t nrows=tmeanROC->GetNrows();
 822     for (UInt_t irow=0;irow<nrows;++irow){
 823       UInt_t npads=tmeanROC->GetNPads(irow);
 824       for (UInt_t ipad=0;ipad<npads;++ipad){
 825         Int_t outlier=0,masked=0;
 826         Float_t q=qmeanROC->GetValue(irow,ipad);
 827         Float_t t=tmeanROC->GetValue(irow,ipad);
 828         //masked channels are outliers
 829         if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
 830         //edge pads are outliers
 831         if (ipad==0||ipad==npads-1) masked=1;
 832         //channels with too large charge or timing deviation from the meadian are outliers
 833 //         if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
 834         if (q<fPulQminLimit && !masked) outlier=1;
 835         //check for nan
 836         if ( !(q<10000000) || !(t<10000000)) outlier=1;
 837         outROC->SetValue(irow,ipad,outlier+masked);
 838         fNpulserOutliers+=outlier;
 839       }
 840     }
 841   }
 842 }
 843 //_____________________________________________________________________________________
 844 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
 845 {
 846   //
 847   // Create pad time0 object from pulser and/or CE data, depending on the selected model
 848   // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
 849   // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
 850   // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
 851   //
 852   // In case model 2 is invoked - gy arival time gradient is also returned
 853   //
 854   gyA=0;
 855   gyC=0;
 856   AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
 857   // decide between different models
 858   if (model==0||model==1){
 859     TVectorD vMean;
 860     if (model==1) ProcessPulser(vMean);
 861     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 862       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 863       if (!rocPulTmean) continue;
 864       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 865       AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
 866       Float_t mean=rocPulTmean->GetMean(rocOut);
 867       //treat case where a whole partition is masked
 868       if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
 869       if (model==1) {
 870         Int_t type=isec/18;
 871         mean=vMean[type];
 872       }
 873       UInt_t nrows=rocTime0->GetNrows();
 874       for (UInt_t irow=0;irow<nrows;++irow){
 875         UInt_t npads=rocTime0->GetNPads(irow);
 876         for (UInt_t ipad=0;ipad<npads;++ipad){
 877           Float_t time=rocPulTmean->GetValue(irow,ipad);
 878           //in case of an outlier pad use the mean of the altro values.
 879           //This should be the most precise guess in that case.
 880           if (rocOut->GetValue(irow,ipad)) {
 881             time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
 882             if ( TMath::Abs(time)<kAlmost0 ) time=mean;
 883           }
 884           Float_t val=time-mean;
 885           rocTime0->SetValue(irow,ipad,val);
 886         }
 887       }
 888     }
 889   } else if (model==2){
 890     Double_t pgya,pgyc,pchi2a,pchi2c;
 891     AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
 892     fCETmean->Add(padPulser,-1.);
 893     TVectorD vA,vC;
 894     AliTPCCalPad outCE("outCE","outCE");
 895     Int_t nOut;
 896     ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
 897     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
 898 //     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
 899     if (!padFit) { delete padPulser; return 0;}
 900     gyA=vA[2];
 901     gyC=vC[2];
 902     fCETmean->Add(padPulser,1.);
 903     padTime0->Add(fCETmean);
 904     padTime0->Add(padFit,-1);
 905     delete padPulser;
 906     TVectorD vFitROC;
 907     TMatrixD mFitROC;
 908     Float_t chi2;
 909     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 910       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 911       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 912       AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
 913       AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
 914       rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
 915       AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
 916       Float_t mean=rocPulTmean->GetMean(rocOutPul);
 917       if ( TMath::Abs(mean)<kAlmost0 ) mean=rocPulTmean->GetMean();
 918       UInt_t nrows=rocTime0->GetNrows();
 919       for (UInt_t irow=0;irow<nrows;++irow){
 920         UInt_t npads=rocTime0->GetNPads(irow);
 921         for (UInt_t ipad=0;ipad<npads;++ipad){
 922           Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
 923           if (rocOutCE->GetValue(irow,ipad)){
 924             Float_t valOut=rocCEfit->GetValue(irow,ipad);
 925             if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
 926             rocTime0->SetValue(irow,ipad,valOut);
 927           }
 928         }
 929       }
 930       delete rocCEfit;
 931     }
 932     delete padFit;
 933   }
 934   Double_t median = padTime0->GetMedian();
 935   padTime0->Add(-median);  // normalize to median
 936   return padTime0;
 937 }
 938 //_____________________________________________________________________________________
 939 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
 940 {
 941   //
 942   // GetMeanAlto information
 943   //
 944   if (roc==0) return 0.;
 945   const Int_t sector=roc->GetSector();
 946   AliTPCROC *tpcRoc=AliTPCROC::Instance();
 947   const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
 948   Float_t mean=0;
 949   Int_t   n=0;
 950
 951   //loop over a small range around the requested pad (+-10 rows/pads)
 952   for (Int_t irow=row-10;irow<row+10;++irow){
 953     if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
 954     for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
 955       if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
 956       const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
 957       if (altroRoc!=altroCurr) continue;
 958       if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
 959       Float_t val=roc->GetValue(irow,ipad);
 960       mean+=val;
 961       ++n;
 962     }
 963   }
 964   if (n>0) mean/=n;
 965   return mean;
 966 }
 967 //_____________________________________________________________________________________
 968 void AliTPCcalibDButil::SetRefFile(const char* filename)
 969 {
 970   //
 971   // load cal pad objects form the reference file
 972   //
 973   TDirectory *currDir=gDirectory;
 974   TFile f(filename);
 975   fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
 976   fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
 977   //pulser data
 978   fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
 979   fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
 980   fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
 981   //CE data
 982   fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
 983   fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
 984   fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
 985   //Altro data
 986 //   fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
 987 //   fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
 988 //   fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
 989 //   fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
 990   fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
 991   f.Close();
 992   currDir->cd();
 993 }
 994 //_____________________________________________________________________________________
 995 void AliTPCcalibDButil::UpdateRefDataFromOCDB()
 996 {
 997   //
 998   // set reference data from OCDB Reference map
 999   //
1000   if (!fRefMap) {
1001     AliWarning("Referenc map not set!");
1002     return;
1003   }
1004
1005   TString cdbPath;
1006   AliCDBEntry* entry = 0x0;
1007   Bool_t hasAnyChanged=kFALSE;
1008
1009   //pedestals
1010   cdbPath="TPC/Calib/Pedestals";
1011   if (HasRefChanged(cdbPath.Data())){
1012     hasAnyChanged=kTRUE;
1013     //delete old entries
1014     if (fRefPedestals) delete fRefPedestals;
1015     if (fRefPedestalMasked) delete fRefPedestalMasked;
1016     fRefPedestals=fRefPedestalMasked=0x0;
1017     //get new entries
1018     entry=GetRefEntry(cdbPath.Data());
1019     if (entry){
1020       entry->SetOwner(kTRUE);
1021       fRefPedestals=GetRefCalPad(entry);
1022       delete entry;
1023       fRefPedestalMasked=GetAltroMasked(cdbPath, "MaskedPedestals");
1024     }
1025   }
1026
1027   //noise
1028   cdbPath="TPC/Calib/PadNoise";
1029   if (HasRefChanged(cdbPath.Data())){
1030     hasAnyChanged=kTRUE;
1031     //delete old entry
1032     if (fRefPadNoise) delete fRefPadNoise;
1033     fRefPadNoise=0x0;
1034     //get new entry
1035     entry=GetRefEntry(cdbPath.Data());
1036     if (entry){
1037       entry->SetOwner(kTRUE);
1038       fRefPadNoise=GetRefCalPad(entry);
1039       delete entry;
1040     }
1041   }
1042
1043   //pulser
1044   cdbPath="TPC/Calib/Pulser";
1045   if (HasRefChanged(cdbPath.Data())){
1046     hasAnyChanged=kTRUE;
1047     //delete old entries
1048     if (fRefPulserTmean) delete fRefPulserTmean;
1049     if (fRefPulserTrms) delete fRefPulserTrms;
1050     if (fRefPulserQmean) delete fRefPulserQmean;
1051     if (fRefPulserMasked) delete fRefPulserMasked;
1052     fRefPulserTmean=fRefPulserTrms=fRefPulserQmean=fRefPulserMasked=0x0;
1053     //get new entries
1054     entry=GetRefEntry(cdbPath.Data());
1055     if (entry){
1056       entry->SetOwner(kTRUE);
1057       fRefPulserTmean=GetRefCalPad(entry,"PulserTmean");
1058       fRefPulserTrms=GetRefCalPad(entry,"PulserTrms");
1059       fRefPulserQmean=GetRefCalPad(entry,"PulserQmean");
1060       delete entry;
1061       fRefPulserMasked=GetAltroMasked(cdbPath, "MaskedPulser");
1062     }
1063   }
1064
1065   //ce
1066   cdbPath="TPC/Calib/CE";
1067   if (HasRefChanged(cdbPath.Data())){
1068     hasAnyChanged=kTRUE;
1069     //delete old entries
1070     if (fRefCETmean) delete fRefCETmean;
1071     if (fRefCETrms) delete fRefCETrms;
1072     if (fRefCEQmean) delete fRefCEQmean;
1073     if (fRefCEMasked) delete fRefCEMasked;
1074     fRefCETmean=fRefCETrms=fRefCEQmean=fRefCEMasked=0x0;
1075     //get new entries
1076     entry=GetRefEntry(cdbPath.Data());
1077     if (entry){
1078       entry->SetOwner(kTRUE);
1079       fRefCETmean=GetRefCalPad(entry,"CETmean");
1080       fRefCETrms=GetRefCalPad(entry,"CETrms");
1081       fRefCEQmean=GetRefCalPad(entry,"CEQmean");
1082       delete entry;
1083       fRefCEMasked=GetAltroMasked(cdbPath, "MaskedCE");
1084     }
1085   }
1086
1087   //altro data
1088   cdbPath="TPC/Calib/AltroConfig";
1089   if (HasRefChanged(cdbPath.Data())){
1090     hasAnyChanged=kTRUE;
1091     //delete old entries
1092     if (fRefALTROFPED) delete fRefALTROFPED;
1093     if (fRefALTROZsThr) delete fRefALTROZsThr;
1094     if (fRefALTROAcqStart) delete fRefALTROAcqStart;
1095     if (fRefALTROAcqStop) delete fRefALTROAcqStop;
1096     if (fRefALTROMasked) delete fRefALTROMasked;
1097     fRefALTROFPED=fRefALTROZsThr=fRefALTROAcqStart=fRefALTROAcqStop=fRefALTROMasked=0x0;
1098     //get new entries
1099     entry=GetRefEntry(cdbPath.Data());
1100     if (entry){
1101       entry->SetOwner(kTRUE);
1102       fRefALTROFPED=GetRefCalPad(entry,"FPED");
1103       fRefALTROZsThr=GetRefCalPad(entry,"ZsThr");
1104       fRefALTROAcqStart=GetRefCalPad(entry,"AcqStart");
1105       fRefALTROAcqStop=GetRefCalPad(entry,"AcqStop");
1106       fRefALTROMasked=GetRefCalPad(entry,"Masked");
1107       delete entry;
1108     }
1109   }
1110
1111   //raw data
1112   /*
1113   cdbPath="TPC/Calib/Raw";
1114   if (HasRefChanged(cdbPath.Data())){
1115     hasAnyChanged=kTRUE;
1116     //delete old entry
1117     if (fRefCalibRaw) delete fRefCalibRaw;
1118     //get new entry
1119     entry=GetRefEntry(cdbPath.Data());
1120     if (entry){
1121       entry->SetOwner(kTRUE);
1122       TObjArray *arr=(TObjArray*)entry->GetObject();
1123       if (!arr){
1124         AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1125       } else {
1126         fRefCalibRaw=(AliTPCCalibRaw*)arr->At(0)->Clone();
1127       }
1128     }
1129   }
1130   */
1131
1132   //data qa
1133   cdbPath="TPC/Calib/QA";
1134   if (HasRefChanged(cdbPath.Data())){
1135     hasAnyChanged=kTRUE;
1136     //delete old entry
1137     if (fRefDataQA) delete fRefDataQA;
1138     //get new entry
1139     entry=GetRefEntry(cdbPath.Data());
1140     if (entry){
1141       entry->SetOwner(kTRUE);
1142       fDataQA=dynamic_cast<AliTPCdataQA*>(entry->GetObject());
1143       if (!fDataQA){
1144         AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1145       } else {
1146         fRefDataQA=(AliTPCdataQA*)fDataQA->Clone();
1147       }
1148       delete entry;
1149     }
1150   }
1151
1152
1153 //update current reference maps
1154   if (hasAnyChanged){
1155     if (fCurrentRefMap) delete fCurrentRefMap;
1156     fCurrentRefMap=(TMap*)fRefMap->Clone();
1157   }
1158 }
1159 //_____________________________________________________________________________________
1160 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry, const char* objName)
1161 {
1162   //
1163   // TObjArray object type case
1164   // find 'objName' in 'arr' cast is to a calPad and store it in 'pad'
1165   //
1166   AliTPCCalPad *pad=0x0;
1167   TObjArray *arr=(TObjArray*)entry->GetObject();
1168   if (!arr){
1169     AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1170     return pad;
1171   }
1172   pad=(AliTPCCalPad*)arr->FindObject(objName);
1173   if (!pad) {
1174     AliError(Form("Could not get '%s' from TObjArray in entry '%s'\nPlease check!!!",objName,entry->GetId().GetPath().Data()));
1175     return pad;
1176   }
1177   return (AliTPCCalPad*)pad->Clone();
1178 }
1179 //_____________________________________________________________________________________
1180 AliTPCCalPad* AliTPCcalibDButil::GetRefCalPad(AliCDBEntry *entry)
1181 {
1182   //
1183   // AliTPCCalPad object type case
1184   // cast object to a calPad and store it in 'pad'
1185   //
1186   AliTPCCalPad *pad=(AliTPCCalPad*)entry->GetObject();
1187   if (!pad) {
1188     AliError(Form("Could not get object from entry '%s'\nPlease check!!!",entry->GetId().GetPath().Data()));
1189     return 0x0;
1190   }
1191   pad=(AliTPCCalPad*)pad->Clone();
1192   return pad;
1193 }
1194 //_____________________________________________________________________________________
1195 AliTPCCalPad* AliTPCcalibDButil::GetAltroMasked(const char* cdbPath, const char* name)
1196 {
1197   //
1198   // set altro masked channel map for 'cdbPath'
1199   //
1200   AliTPCCalPad* pad=0x0;
1201   const Int_t run=GetReferenceRun(cdbPath);
1202   if (run<0) {
1203     AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1204     return pad;
1205   }
1206   AliCDBEntry *entry=AliCDBManager::Instance()->Get("TPC/Calib/AltroConfig", run);
1207   if (!entry) {
1208     AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1209     return pad;
1210   }
1211   pad=GetRefCalPad(entry,"Masked");
1212   if (pad) pad->SetNameTitle(name,name);
1213   entry->SetOwner(kTRUE);
1214   delete entry;
1215   return pad;
1216 }
1217 //_____________________________________________________________________________________
1218 void AliTPCcalibDButil::SetReferenceRun(Int_t run){
1219   //
1220   // Get Reference map
1221   //
1222   if (run<0) run=fCalibDB->GetRun();
1223   TString cdbPath="TPC/Calib/Ref";
1224   AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath.Data(), run);
1225   if (!entry) {
1226     AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath.Data()));
1227     fRefMap=0;
1228     return;
1229   }
1230   entry->SetOwner(kTRUE);
1231   fRefMap=(TMap*)(entry->GetObject());
1232   AliCDBId &id=entry->GetId();
1233   fRefValidity.Form("%d_%d_v%d_s%d",id.GetFirstRun(),id.GetLastRun(),id.GetVersion(),id.GetSubVersion());
1234 }
1235 //_____________________________________________________________________________________
1236 Bool_t AliTPCcalibDButil::HasRefChanged(const char *cdbPath)
1237 {
1238   //
1239   // check whether a reference cdb entry has changed
1240   //
1241   if (!fCurrentRefMap) return kTRUE;
1242   if (GetReferenceRun(cdbPath)!=GetCurrentReferenceRun(cdbPath)) return kTRUE;
1243   return kFALSE;
1244 }
1245 //_____________________________________________________________________________________
1246 AliCDBEntry* AliTPCcalibDButil::GetRefEntry(const char* cdbPath)
1247 {
1248   //
1249   // get the reference AliCDBEntry for 'cdbPath'
1250   //
1251   const Int_t run=GetReferenceRun(cdbPath);
1252   if (run<0) {
1253     AliError(Form("Could not get reference run number for object '%s'\nPlease check availability!!!",cdbPath));
1254     return 0;
1255   }
1256   AliCDBEntry *entry=AliCDBManager::Instance()->Get(cdbPath, run);
1257   if (!entry) {
1258     AliError(Form("Could not get reference object '%s'\nPlease check availability!!!",cdbPath));
1259     return 0;
1260   }
1261   return entry;
1262 }
1263 //_____________________________________________________________________________________
1264 const Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type){
1265   //
1266   // Get reference run number for the specified OCDB path
1267   //
1268   if (!fCurrentRefMap) return -2;
1269   TObjString *str=dynamic_cast<TObjString*>(fCurrentRefMap->GetValue(type));
1270   if (!str) return -2;
1271   return (const Int_t)str->GetString().Atoi();
1272 }
1273 //_____________________________________________________________________________________
1274 const Int_t AliTPCcalibDButil::GetReferenceRun(const char* type) const{
1275   //
1276   // Get reference run number for the specified OCDB path
1277   //
1278   if (!fRefMap) return -1;
1279   TObjString *str=dynamic_cast<TObjString*>(fRefMap->GetValue(type));
1280   if (!str) return -1;
1281   return (const Int_t)str->GetString().Atoi();
1282 }
1283 //_____________________________________________________________________________________
1284 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad * const ceOut, Float_t minSignal, Float_t cutTrmsMin,  Float_t cutTrmsMax, Float_t cutMaxDistT){
1285   //
1286   // Author:  marian.ivanov@cern.ch
1287   //
1288   // Create outlier map for CE study
1289   // Parameters:
1290   //  Return value - outlyer map
1291   //  noutlyersCE  - number of outlyers
1292   //  minSignal    - minimal total Q signal
1293   //  cutRMSMin    - minimal width of the signal in respect to the median
1294   //  cutRMSMax    - maximal width of the signal in respect to the median
1295   //  cutMaxDistT  - maximal deviation from time median per chamber
1296   //
1297   // Outlyers criteria:
1298   // 0. Exclude masked pads
1299   // 1. Exclude first two rows in IROC and last two rows in OROC
1300   // 2. Exclude edge pads
1301   // 3. Exclude channels with too large variations
1302   // 4. Exclude pads with too small signal
1303   // 5. Exclude signal with outlyers RMS
1304   // 6. Exclude channels to far from the chamber median
1305   noutliersCE=0;
1306   //create outlier map
1307   AliTPCCalPad *out=ceOut;
1308   if (!out)     out= new AliTPCCalPad("outCE","outCE");
1309   AliTPCCalROC *rocMasked=0x0;
1310   if (!fCETmean) return 0;
1311   if (!fCETrms) return 0;
1312   if (!fCEQmean) return 0;
1313   //
1314   //loop over all channels
1315   //
1316   Double_t rmsMedian         = fCETrms->GetMedian();
1317   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1318     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
1319     if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1320     AliTPCCalROC *rocOut       = out->GetCalROC(iroc);
1321     AliTPCCalROC *rocCEQ       = fCEQmean->GetCalROC(iroc);
1322     AliTPCCalROC *rocCETrms    = fCETrms->GetCalROC(iroc);
1323     Double_t trocMedian        = rocData->GetMedian();
1324     //
1325     if (!rocData) {
1326       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1327       rocOut->Add(1.);
1328       continue;
1329     }
1330     //
1331     //select outliers
1332     UInt_t nrows=rocData->GetNrows();
1333     for (UInt_t irow=0;irow<nrows;++irow){
1334       UInt_t npads=rocData->GetNPads(irow);
1335       for (UInt_t ipad=0;ipad<npads;++ipad){
1336         rocOut->SetValue(irow,ipad,0);
1337         Float_t valTmean=rocData->GetValue(irow,ipad);
1338         Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1339         Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1340         //0. exclude masked pads
1341         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1342           rocOut->SetValue(irow,ipad,1);
1343           continue;
1344         }
1345         //1. exclude first two rows in IROC and last two rows in OROC
1346         if (iroc<36){
1347           if (irow<2) rocOut->SetValue(irow,ipad,1);
1348         } else {
1349           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1350         }
1351         //2. exclude edge pads
1352         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1353         //exclude values that are exactly 0
1354         if ( TMath::Abs(valTmean)<kAlmost0) {
1355           rocOut->SetValue(irow,ipad,1);
1356           ++noutliersCE;
1357         }
1358         //3.  exclude channels with too large variations
1359         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1360           rocOut->SetValue(irow,ipad,1);
1361           ++noutliersCE;
1362         }
1363         //
1364         //4.  exclude channels with too small signal
1365         if (valQmean<minSignal) {
1366           rocOut->SetValue(irow,ipad,1);
1367           ++noutliersCE;
1368         }
1369         //
1370         //5. exclude channels with too small rms
1371         if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1372           rocOut->SetValue(irow,ipad,1);
1373           ++noutliersCE;
1374         }
1375         //
1376         //6. exclude channels to far from the chamber median
1377         if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1378           rocOut->SetValue(irow,ipad,1);
1379           ++noutliersCE;
1380         }
1381       }
1382     }
1383   }
1384   //
1385   return out;
1386 }
1387
1388
1389 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1390   //
1391   // Author: marian.ivanov@cern.ch
1392   //
1393   // Create outlier map for Pulser
1394   // Parameters:
1395   //  Return value     - outlyer map
1396   //  noutlyersPulser  - number of outlyers
1397   //  cutTime          - absolute cut - distance to the median of chamber
1398   //  cutnRMSQ         - nsigma cut from median  q distribution per chamber
1399   //  cutnRMSrms       - nsigma cut from median  rms distribution
1400   // Outlyers criteria:
1401   // 0. Exclude masked pads
1402   // 1. Exclude time outlyers (default 3 time bins)
1403   // 2. Exclude q outlyers    (default 5 sigma)
1404   // 3. Exclude rms outlyers  (default 5 sigma)
1405   noutliersPulser=0;
1406   AliTPCCalPad *out=pulserOut;
1407   if (!out)     out= new AliTPCCalPad("outPulser","outPulser");
1408   AliTPCCalROC *rocMasked=0x0;
1409   if (!fPulserTmean) return 0;
1410   if (!fPulserTrms) return 0;
1411   if (!fPulserQmean) return 0;
1412   //
1413   //loop over all channels
1414   //
1415   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1416     if (fALTROMasked)   rocMasked= fALTROMasked->GetCalROC(iroc);
1417     AliTPCCalROC *rocData       = fPulserTmean->GetCalROC(iroc);
1418     AliTPCCalROC *rocOut        = out->GetCalROC(iroc);
1419     AliTPCCalROC *rocPulserQ    = fPulserQmean->GetCalROC(iroc);
1420     AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1421     //
1422     Double_t rocMedianT         = rocData->GetMedian();
1423     Double_t rocMedianQ         = rocPulserQ->GetMedian();
1424     Double_t rocRMSQ            = rocPulserQ->GetRMS();
1425     Double_t rocMedianTrms      = rocPulserTrms->GetMedian();
1426     Double_t rocRMSTrms         = rocPulserTrms->GetRMS();
1427     for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1428       rocOut->SetValue(ichannel,0);
1429       Float_t valTmean=rocData->GetValue(ichannel);
1430       Float_t valQmean=rocPulserQ->GetValue(ichannel);
1431       Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1432       Int_t isOut=0;
1433       if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1434       if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1435       if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1436       rocOut->SetValue(ichannel,isOut);
1437       if (isOut) noutliersPulser++;
1438     }
1439   }
1440   return out;
1441 }
1442
1443
1444 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1445   //
1446   // Author : Marian Ivanov
1447   // Create pad time0 correction map using information from the CE and from pulser
1448   //
1449   //
1450   // Return PadTime0 to be used for time0 relative alignment
1451   // if dump file specified intermediat results are dumped to the fiel and can be visualized
1452   // using $ALICE_ROOT/TPC/script/gui application
1453   //
1454   // fitResultsA - fitParameters A side
1455   // fitResultsC - fitParameters C side
1456   // chi2A       - chi2/ndf for A side (assuming error 1 time bin)
1457   // chi2C       - chi2/ndf for C side (assuming error 1 time bin)
1458   //
1459   //
1460   // Algorithm:
1461   // 1. Find outlier map for CE
1462   // 2. Find outlier map for Pulser
1463   // 3. Replace outlier by median at given sector  (median without outliers)
1464   // 4. Substract from the CE data pulser
1465   // 5. Fit the CE with formula
1466   //    5.1) (IROC-OROC) offset
1467   //    5.2) gx
1468   //    5.3) gy
1469   //    5.4) (lx-xmid)
1470   //    5.5) (IROC-OROC)*(lx-xmid)
1471   //    5.6) (ly/lx)^2
1472   // 6. Substract gy fit dependence from the CE data
1473   // 7. Add pulser back to CE data
1474   // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1475   // 9. return CE data
1476   //
1477   // Time0 <= padCE = padCEin  -padCEfitGy  - if not outlier
1478   // Time0 <= padCE = padFitAll-padCEfitGy  - if outlier
1479
1480   // fit formula
1481   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)";
1482   // output for fit formula
1483   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)";
1484   // gy part of formula
1485   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)";
1486   //
1487   //
1488   if (!fCETmean) return 0;
1489   Double_t pgya,pgyc,pchi2a,pchi2c;
1490   AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1491   AliTPCCalPad * padCEOut     = CreateCEOutlyerMap(nOut);
1492
1493   AliTPCCalPad * padPulser    = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1494   AliTPCCalPad * padCE        = new AliTPCCalPad(*fCETmean);
1495   AliTPCCalPad * padCEIn      = new AliTPCCalPad(*fCETmean);
1496   AliTPCCalPad * padOut       = new AliTPCCalPad("padOut","padOut");
1497   padPulser->SetName("padPulser");
1498   padPulserOut->SetName("padPulserOut");
1499   padCE->SetName("padCE");
1500   padCEIn->SetName("padCEIn");
1501   padCEOut->SetName("padCEOut");
1502   padOut->SetName("padOut");
1503
1504   //
1505   // make combined outlyers map
1506   // and replace outlyers in maps with median for chamber
1507   //
1508   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1509     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1510     AliTPCCalROC * rocPulser    = padPulser->GetCalROC(iroc);
1511     AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1512     AliTPCCalROC * rocCEOut     = padCEOut->GetCalROC(iroc);
1513     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1514     Double_t ceMedian           = rocCE->GetMedian(rocCEOut);
1515     Double_t pulserMedian       = rocPulser->GetMedian(rocCEOut);
1516     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1517       if (rocPulserOut->GetValue(ichannel)>0) {
1518         rocPulser->SetValue(ichannel,pulserMedian);
1519         rocOut->SetValue(ichannel,1);
1520       }
1521       if (rocCEOut->GetValue(ichannel)>0) {
1522         rocCE->SetValue(ichannel,ceMedian);
1523         rocOut->SetValue(ichannel,1);
1524       }
1525     }
1526   }
1527   //
1528   // remove pulser time 0
1529   //
1530   padCE->Add(padPulser,-1);
1531   //
1532   // Make fits
1533   //
1534   TMatrixD dummy;
1535   Float_t chi2Af,chi2Cf;
1536   padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1537   chi2A=chi2Af;
1538   chi2C=chi2Cf;
1539   //
1540   AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1541   padCEFitGY->SetName("padCEFitGy");
1542   //
1543   AliTPCCalPad *padCEFit  =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1544   padCEFit->SetName("padCEFit");
1545   //
1546   AliTPCCalPad* padCEDiff  = new AliTPCCalPad(*padCE);
1547   padCEDiff->SetName("padCEDiff");
1548   padCEDiff->Add(padCEFit,-1.);
1549   //
1550   //
1551   padCE->Add(padCEFitGY,-1.);
1552
1553   padCE->Add(padPulser,1.);
1554   Double_t padmedian = padCE->GetMedian();
1555   padCE->Add(-padmedian);  // normalize to median
1556   //
1557   // Replace outliers by fit value - median of diff per given chamber -GY fit
1558   //
1559   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1560     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1561     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1562     AliTPCCalROC * rocCEFit     = padCEFit->GetCalROC(iroc);
1563     AliTPCCalROC * rocCEFitGY   = padCEFitGY->GetCalROC(iroc);
1564     AliTPCCalROC * rocCEDiff    = padCEDiff->GetCalROC(iroc);
1565     //
1566     Double_t diffMedian         = rocCEDiff->GetMedian(rocOut);
1567     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1568       if (rocOut->GetValue(ichannel)==0) continue;
1569       Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1570       rocCE->SetValue(ichannel,value);
1571     }
1572   }
1573   //
1574   //
1575   if (dumpfile){
1576     //dump to the file - result can be visualized
1577     AliTPCPreprocessorOnline preprocesor;
1578     preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1579     preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1580     preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1581     preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1582     //
1583     preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1584     preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1585     //
1586     preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1587     preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1588     preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1589     preprocesor.DumpToFile(dumpfile);
1590   }
1591   delete padPulser;
1592   delete padPulserOut;
1593   delete padCEIn;
1594   delete padCEOut;
1595   delete padOut;
1596   delete padCEDiff;
1597   delete padCEFitGY;
1598   return padCE;
1599 }
1600
1601
1602
1603
1604
1605 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1606   //
1607   // find the closest point to xref  in x  direction
1608   // return dx and value
1609   Int_t index=0;
1610   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1611   if (index<0) index=0;
1612   if (index>=graph->GetN()-1) index=graph->GetN()-2;
1613   if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1614   dx = xref-graph->GetX()[index];
1615   y  = graph->GetY()[index];
1616   return index;
1617 }
1618
1619
1620 Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1621   //
1622   // Get the correction of the trigger offset
1623   // combining information from the laser track calibration
1624   // and from cosmic calibration
1625   //
1626   // run       - run number
1627   // timeStamp - tim stamp in seconds
1628   // deltaT    - integration period to calculate offset
1629   // deltaTLaser -max validity of laser data
1630   // valType   - 0 - median, 1- mean
1631   //
1632   // Integration vaues are just recomendation - if not possible to get points
1633   // automatically increase the validity by factor 2
1634   // (recursive algorithm until one month of data taking)
1635   //
1636   //
1637   const Float_t kLaserCut=0.0005;
1638   const Int_t   kMaxPeriod=3600*24*30*3; // 3 month max
1639   const Int_t   kMinPoints=20;
1640   //
1641   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1642   if (!array) {
1643     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1644   }
1645   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1646   if (!array) return 0;
1647   //
1648   TGraphErrors *laserA[3]={0,0,0};
1649   TGraphErrors *laserC[3]={0,0,0};
1650   TGraphErrors *cosmicAll=0;
1651   laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1652   laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1653   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1654   //
1655   //
1656   if (!cosmicAll) return 0;
1657   Int_t nmeasC=cosmicAll->GetN();
1658   Float_t *tdelta = new Float_t[nmeasC];
1659   Int_t nused=0;
1660   for (Int_t i=0;i<nmeasC;i++){
1661     if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1662     Float_t ccosmic=cosmicAll->GetY()[i];
1663     Double_t yA=0,yC=0,dA=0,dC=0;
1664     if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1665     if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1666     //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1667     //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1668     //
1669     if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1670     Float_t claser=0;
1671     if (TMath::Abs(yA-yC)<kLaserCut) {
1672       claser=(yA-yC)*0.5;
1673     }else{
1674       if (i%2==0)  claser=yA;
1675       if (i%2==1)  claser=yC;
1676     }
1677     tdelta[nused]=ccosmic-claser;
1678     nused++;
1679   }
1680   if (nused<kMinPoints &&deltaT<kMaxPeriod) return  AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1681   Double_t median = TMath::Median(nused,tdelta);
1682   Double_t mean  = TMath::Mean(nused,tdelta);
1683   delete tdelta;
1684   return (valType==0) ? median:mean;
1685 }
1686
1687 Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1688   //
1689   // Get the correction of the drift velocity
1690   // combining information from the laser track calibration
1691   // and from cosmic calibration
1692   //
1693   // dist      - return value - distance to closest point in graph
1694   // run       - run number
1695   // timeStamp - tim stamp in seconds
1696   // deltaT    - integration period to calculate time0 offset
1697   // deltaTLaser -max validity of laser data
1698   // valType   - 0 - median, 1- mean
1699   //
1700   // Integration vaues are just recomendation - if not possible to get points
1701   // automatically increase the validity by factor 2
1702   // (recursive algorithm until one month of data taking)
1703   //
1704   //
1705   //
1706   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1707   if (!array) {
1708     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1709   }
1710   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1711   if (!array) return 0;
1712   TGraphErrors *cosmicAll=0;
1713   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1714   if (!cosmicAll) return 0;
1715   Double_t grY=0;
1716   AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1717
1718   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1719   Double_t vcosmic=  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1720   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1721   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
1722   return  vcosmic-t0;
1723
1724   /*
1725     Example usage:
1726
1727     Int_t run=89000
1728     TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1729     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1730     laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1731     //
1732     Double_t *yvd= new Double_t[cosmicAll->GetN()];
1733     Double_t *yt0= new Double_t[cosmicAll->GetN()];
1734     for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1735     for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1736
1737     TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1738     TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1739
1740   */
1741
1742 }
1743
1744 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1745 {
1746   //
1747   // Create a default name for the gui file
1748   //
1749
1750   return Form("guiRefTreeRun%s.root",GetRefValidity());
1751 }
1752
1753 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1754 {
1755   //
1756   // Create a gui reference tree
1757   // if dirname and filename are empty default values will be used
1758   // this is the recommended way of using this function
1759   // it allows to check whether a file with the given run validity alredy exists
1760   //
1761   if (!AliCDBManager::Instance()->GetDefaultStorage()){
1762     AliError("Default Storage not set. Cannot create reference calibration Tree!");
1763     return kFALSE;
1764   }
1765
1766   TString file=filename;
1767   if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1768
1769   AliTPCPreprocessorOnline prep;
1770   //noise and pedestals
1771   if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1772   if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1773   if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1774   //pulser data
1775   if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1776   if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1777   if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1778   if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1779   //CE data
1780   if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1781   if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1782   if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1783   if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1784   //Altro data
1785   if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1786   if (fRefALTROZsThr    ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr    )));
1787   if (fRefALTROFPED     ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED     )));
1788   if (fRefALTROAcqStop  ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop  )));
1789   if (fRefALTROMasked   ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked   )));
1790   //QA
1791   AliTPCdataQA *dataQA=fRefDataQA;
1792   if (dataQA) {
1793     if (dataQA->GetNLocalMaxima())
1794       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1795     if (dataQA->GetMaxCharge())
1796       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1797     if (dataQA->GetMeanCharge())
1798       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1799     if (dataQA->GetNoThreshold())
1800       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1801     if (dataQA->GetNTimeBins())
1802       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1803     if (dataQA->GetNPads())
1804       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1805     if (dataQA->GetTimePosition())
1806       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1807   }
1808   prep.DumpToFile(file.Data());
1809   return kTRUE;
1810 }
1811
1812 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1813   //
1814   // Get the correction of the drift velocity using the laser tracks calbration
1815   //
1816   // run       - run number
1817   // timeStamp - tim stamp in seconds
1818   // deltaT    - integration period to calculate time0 offset
1819   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1820   // Note in case no data form both A and C side - the value from active side used
1821   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1822   TGraphErrors *grlaserA=0;
1823   TGraphErrors *grlaserC=0;
1824   Double_t vlaserA=0, vlaserC=0;
1825   if (!array) return 0;
1826   grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1827   grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1828   Double_t deltaY;
1829   if (grlaserA) {
1830     AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1831     if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
1832     else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1833   }
1834   if (grlaserC) {
1835     AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1836     if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
1837     else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1838   }
1839   if (side==0) return vlaserA;
1840   if (side==1) return vlaserC;
1841   Double_t mdrift=(vlaserA+vlaserC)*0.5;
1842   if (!grlaserA) return vlaserC;
1843   if (!grlaserC) return vlaserA;
1844   return mdrift;
1845 }
1846
1847
1848 Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1849   //
1850   // Get the correction of the drift velocity using the CE laser data
1851   // combining information from the CE,  laser track calibration
1852   // and P/T calibration
1853   //
1854   // run       - run number
1855   // timeStamp - tim stamp in seconds
1856   // deltaT    - integration period to calculate time0 offset
1857   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1858   TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
1859   if (!arrT) return 0;
1860   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
1861   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1862   AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1863   //
1864   //
1865   Double_t corrPTA = 0, corrPTC=0;
1866   Double_t ltime0A = 0, ltime0C=0;
1867   Double_t gry=0;
1868   Double_t corrA=0, corrC=0;
1869   Double_t timeA=0, timeC=0;
1870   TGraph *graphA = (TGraph*)arrT->At(72);
1871   TGraph *graphC = (TGraph*)arrT->At(73);
1872   if (!graphA && !graphC) return 0.;
1873   if (graphA &&graphA->GetN()>0) {
1874     AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1875     timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1876     Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1877     ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1878     if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
1879     corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1880     corrA-=corrPTA;
1881   }
1882   if (graphC&&graphC->GetN()>0){
1883     AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1884     timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1885     Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1886     ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1887     if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
1888     corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1889     corrC-=corrPTC;
1890   }
1891
1892   if (side ==0 ) return corrA;
1893   if (side ==1 ) return corrC;
1894   Double_t corrM= (corrA+corrC)*0.5;
1895   if (!graphA) corrM=corrC;
1896   if (!graphC) corrM=corrA;
1897   return corrM;
1898 }
1899
1900 Double_t  AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1901   //
1902   // return drift velocity using the TPC-ITS matchin method
1903   // return also distance to the closest point
1904   //
1905   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1906   TGraphErrors *graph=0;
1907   dist=0;
1908   if (!array) return 0;
1909   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
1910   if (!graph) return 0;
1911   Double_t deltaY;
1912   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1913   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1914   return value;
1915 }
1916
1917 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1918   //
1919   // Get time dependent time 0 (trigger delay in cm) correction
1920   // Arguments:
1921   // timestamp - timestamp
1922   // run       - run number
1923   //
1924   // Notice - Extrapolation outside of calibration range  - using constant function
1925   //
1926   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1927   TGraphErrors *graph=0;
1928   dist=0;
1929   if (!array) return 0;
1930   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
1931   if (!graph) return 0;
1932   Double_t deltaY;
1933   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1934   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1935   return value;
1936 }
1937
1938
1939
1940
1941
1942 Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1943   //
1944   // VERY obscure method - we need something in framework
1945   // Find the TPC runs with temperature OCDB entry
1946   // cache the start and end of the run
1947   //
1948   AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1949   if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1950   if (!storage) return 0;
1951   TString path=storage->GetURI();
1952   TString runsT;
1953   {
1954     TString command;
1955     if (path.Contains("local")){  // find the list if local system
1956       path.ReplaceAll("local://","");
1957       path+="TPC/Calib/Temperature";
1958       command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
1959     }
1960     runsT=gSystem->GetFromPipe(command);
1961   }
1962   TObjArray *arr= runsT.Tokenize("r");
1963   if (!arr) return 0;
1964   //
1965   TArrayI indexes(arr->GetEntries());
1966   TArrayI runs(arr->GetEntries());
1967   Int_t naccept=0;
1968   {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1969       Int_t irunN = atoi(arr->At(irun)->GetName());
1970       if (irunN<startRun) continue;
1971       if (irunN>stopRun) continue;
1972       runs[naccept]=irunN;
1973       naccept++;
1974     }}
1975   fRuns.Set(naccept);
1976   fRunsStart.Set(fRuns.fN);
1977   fRunsStop.Set(fRuns.fN);
1978   TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1979   for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
1980
1981   //
1982   AliCDBEntry * entry = 0;
1983   {for (Int_t irun=0;irun<fRuns.fN; irun++){
1984       entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1985       if (!entry) continue;
1986       AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1987       if (!tmpRun) continue;
1988       fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1989       fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1990       //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1991     }}
1992   return fRuns.fN;
1993 }
1994
1995
1996 Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
1997   //
1998   // binary search - find the run for given time stamp
1999   //
2000   Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2001   Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2002   Int_t cindex  = -1;
2003   for (Int_t index=index0; index<=index1; index++){
2004     if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2005     if (debug) {
2006       printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2007     }
2008   }
2009   if (cindex<0) cindex =(index0+index1)/2;
2010   if (cindex<0) {
2011     return 0;
2012   }
2013   return fRuns[cindex];
2014 }
2015
2016
2017
2018
2019
2020 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2021   //
2022   // filter outlyer measurement
2023   // Only points around median +- sigmaCut filtered
2024   if (!graph) return  0;
2025   Int_t kMinPoints=2;
2026   Int_t npoints0 = graph->GetN();
2027   Int_t npoints=0;
2028   Float_t  rmsY=0;
2029   Double_t *outx=new Double_t[npoints0];
2030   Double_t *outy=new Double_t[npoints0];
2031   //
2032   //
2033   if (npoints0<kMinPoints) return 0;
2034   for (Int_t iter=0; iter<3; iter++){
2035     npoints=0;
2036     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2037       if (graph->GetY()[ipoint]==0) continue;
2038       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2039       outx[npoints]  = graph->GetX()[ipoint];
2040       outy[npoints]  = graph->GetY()[ipoint];
2041       npoints++;
2042     }
2043     if (npoints<=1) break;
2044     medianY  =TMath::Median(npoints,outy);
2045     rmsY   =TMath::RMS(npoints,outy);
2046   }
2047   TGraph *graphOut=0;
2048   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2049   return graphOut;
2050 }
2051
2052
2053 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2054   //
2055   // filter outlyer measurement
2056   // Only points around median +- cut filtered
2057   if (!graph) return  0;
2058   Int_t kMinPoints=2;
2059   Int_t npoints0 = graph->GetN();
2060   Int_t npoints=0;
2061   Float_t  rmsY=0;
2062   Double_t *outx=new Double_t[npoints0];
2063   Double_t *outy=new Double_t[npoints0];
2064   //
2065   //
2066   if (npoints0<kMinPoints) return 0;
2067   for (Int_t iter=0; iter<3; iter++){
2068     npoints=0;
2069     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2070       if (graph->GetY()[ipoint]==0) continue;
2071       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2072       outx[npoints]  = graph->GetX()[ipoint];
2073       outy[npoints]  = graph->GetY()[ipoint];
2074       npoints++;
2075     }
2076     if (npoints<=1) break;
2077     medianY  =TMath::Median(npoints,outy);
2078     rmsY   =TMath::RMS(npoints,outy);
2079   }
2080   TGraph *graphOut=0;
2081   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2082   return graphOut;
2083 }
2084
2085
2086
2087 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2088   //
2089   // filter outlyer measurement
2090   // Only points with normalized errors median +- sigmaCut filtered
2091   //
2092   Int_t kMinPoints=10;
2093   Int_t npoints0 = graph->GetN();
2094   Int_t npoints=0;
2095   Float_t  medianErr=0, rmsErr=0;
2096   Double_t *outx=new Double_t[npoints0];
2097   Double_t *outy=new Double_t[npoints0];
2098   Double_t *erry=new Double_t[npoints0];
2099   Double_t *nerry=new Double_t[npoints0];
2100   Double_t *errx=new Double_t[npoints0];
2101   //
2102   //
2103   if (npoints0<kMinPoints) return 0;
2104   for (Int_t iter=0; iter<3; iter++){
2105     npoints=0;
2106     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2107       nerry[npoints]  = graph->GetErrorY(ipoint);
2108       if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2109       erry[npoints]  = graph->GetErrorY(ipoint);
2110       outx[npoints]  = graph->GetX()[ipoint];
2111       outy[npoints]  = graph->GetY()[ipoint];
2112       errx[npoints]  = graph->GetErrorY(ipoint);
2113       npoints++;
2114     }
2115     if (npoints==0) break;
2116     medianErr=TMath::Median(npoints,erry);
2117     medianY  =TMath::Median(npoints,outy);
2118     rmsErr   =TMath::RMS(npoints,erry);
2119   }
2120   TGraphErrors *graphOut=0;
2121   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2122   delete []outx;
2123   delete []outy;
2124   delete []errx;
2125   delete []erry;
2126   return graphOut;
2127 }
2128
2129
2130 void AliTPCcalibDButil::Sort(TGraph *graph){
2131   //
2132   // sort array - neccessay for approx
2133   //
2134   Int_t npoints = graph->GetN();
2135   Int_t *indexes=new Int_t[npoints];
2136   Double_t *outx=new Double_t[npoints];
2137   Double_t *outy=new Double_t[npoints];
2138   TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2139   for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2140   for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2141   for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2142   for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2143
2144 }
2145 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2146   //
2147   // smmoth graph - mean on the interval
2148   //
2149   Sort(graph);
2150   Int_t npoints = graph->GetN();
2151   Double_t *outy=new Double_t[npoints];
2152
2153   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2154     Double_t lx=graph->GetX()[ipoint];
2155     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2156     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2157     if (index0<0) index0=0;
2158     if (index1>=npoints-1) index1=npoints-1;
2159     if ((index1-index0)>1){
2160       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2161     }else{
2162       outy[ipoint]=graph->GetY()[ipoint];
2163     }
2164   }
2165  //  TLinearFitter  fitter(3,"pol2");
2166 //   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2167 //     Double_t lx=graph->GetX()[ipoint];
2168 //     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2169 //     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2170 //     if (index0<0) index0=0;
2171 //     if (index1>=npoints-1) index1=npoints-1;
2172 //     fitter.ClearPoints();
2173 //     for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2174 //     if ((index1-index0)>1){
2175 //       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2176 //     }else{
2177 //       outy[ipoint]=graph->GetY()[ipoint];
2178 //     }
2179 //   }
2180
2181
2182
2183   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2184     graph->GetY()[ipoint] = outy[ipoint];
2185   }
2186   delete[] outy;
2187 }
2188
2189 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2190   //
2191   // Use constant interpolation outside of range
2192   //
2193   if (!graph) {
2194     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2195     return 0;
2196   }
2197   if (graph->GetN()<1){
2198     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2199     return 0;
2200   }
2201   if (xref<graph->GetX()[0]) return graph->GetY()[0];
2202   if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2203   return graph->Eval( xref);
2204 }
2205
2206 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy,  Double_t sigmaCut){
2207   //
2208   // Filter DCS sensor information
2209   //   ymin     - minimal value
2210   //   ymax     - max value
2211   //   maxdy    - maximal deirivative
2212   //   sigmaCut - cut on values and derivative in terms of RMS distribution
2213   // Return value - accepted fraction
2214   //
2215   // Algorithm:
2216   //
2217   // 0. Calculate median and rms of values in specified range
2218   // 1. Filter out outliers - median+-sigmaCut*rms
2219   //    values replaced by median
2220   //
2221   AliSplineFit * fit    = sensor->GetFit();
2222   if (!fit) return 0.;
2223   Int_t          nknots = fit->GetKnots();
2224   if (nknots==0) {
2225     delete fit;
2226     sensor->SetFit(0);
2227     return 0;
2228   }
2229   //
2230   Double_t *yin0  = new Double_t[nknots];
2231   Double_t *yin1  = new Double_t[nknots];
2232   Int_t naccept=0;
2233
2234   for (Int_t iknot=0; iknot< nknots; iknot++){
2235     if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2236       yin0[naccept]  = fit->GetY0()[iknot];
2237       yin1[naccept]  = fit->GetY1()[iknot];
2238       if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2239       naccept++;
2240     }
2241   }
2242   if (naccept<1) {
2243     delete fit;
2244     sensor->SetFit(0);
2245     return 0.;
2246   }
2247
2248   Double_t medianY0=0, medianY1=0;
2249   Double_t rmsY0   =0, rmsY1=0;
2250   medianY0 = TMath::Median(naccept, yin0);
2251   medianY1 = TMath::Median(naccept, yin1);
2252   rmsY0    = TMath::RMS(naccept, yin0);
2253   rmsY1    = TMath::RMS(naccept, yin1);
2254   naccept=0;
2255   //
2256   // 1. Filter out outliers - median+-sigmaCut*rms
2257   //    values replaced by median
2258   //    if replaced the derivative set to 0
2259   //
2260   for (Int_t iknot=0; iknot< nknots; iknot++){
2261     Bool_t isOK=kTRUE;
2262     if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2263     if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2264     if (nknots<2) fit->GetY1()[iknot]=0;
2265     if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2266     if (!isOK){
2267       fit->GetY0()[iknot]=medianY0;
2268       fit->GetY1()[iknot]=0;
2269     }else{
2270       naccept++;
2271     }
2272   }
2273   delete [] yin0;
2274   delete [] yin1;
2275   return Float_t(naccept)/Float_t(nknots);
2276 }
2277
2278 Float_t  AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2279   //
2280   // Filter temperature array
2281   // tempArray    - array of temperatures         -
2282   // ymin         - minimal accepted temperature  - default 15
2283   // ymax         - maximal accepted temperature  - default 22
2284   // sigmaCut     - values filtered on interval median+-sigmaCut*rms - defaut 5
2285   // return value - fraction of filtered sensors
2286   const Double_t kMaxDy=0.1;
2287   Int_t nsensors=tempArray->NumSensors();
2288   if (nsensors==0) return 0.;
2289   Int_t naccept=0;
2290   for (Int_t isensor=0; isensor<nsensors; isensor++){
2291     AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2292     if (!sensor) continue;
2293     //printf("%d\n",isensor);
2294     FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2295     if (sensor->GetFit()==0){
2296       //delete sensor;
2297       tempArray->RemoveSensorNum(isensor);
2298     }else{
2299       naccept++;
2300     }
2301   }
2302   return Float_t(naccept)/Float_t(nsensors);
2303 }
2304
2305
2306 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2307   //
2308   // Filter CE data
2309   // Input parameters:
2310   //    deltaT   - smoothing window (in seconds)
2311   //    cutAbs   - max distance of the time info to the median (in time bins)
2312   //    cutSigma - max distance (in the RMS)
2313   //    pcstream - optional debug streamer to store original and filtered info
2314   // Hardwired parameters:
2315   //    kMinPoints =10;       // minimal number of points to define the CE
2316   //    kMinSectors=12;       // minimal number of sectors to define sideCE
2317   // Algorithm:
2318   // 0. Filter almost emty graphs (kMinPoints=10)
2319   // 1. calculate median and RMS per side
2320   // 2. Filter graphs - in respect with side medians
2321   //                  - cutAbs and cutDelta used
2322   // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2323   // 4. Calculate mean for A side and C side
2324   //
2325   const Int_t kMinPoints =10;       // minimal number of points to define the CE
2326   const Int_t kMinSectors=12;       // minimal number of sectors to define sideCE
2327   const Int_t kMinTime   =400;     // minimal arrival time of CE
2328   TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2329   Double_t medianY=0;
2330   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
2331   if (!cearray) return;
2332   Double_t tmin=-1;
2333   Double_t tmax=-1;
2334   //
2335   //
2336   AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2337   AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernPressure");
2338   if ( tempMapCE && cavernPressureCE){
2339     //
2340     Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2341     FilterSensor(cavernPressureCE,960,1050,10, 5.);
2342     if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2343     if (isOK)  {
2344       // recalculate P/T correction map for time of the CE
2345       AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2346       driftCalib->SetName("driftPTCE");
2347       driftCalib->SetTitle("driftPTCE");
2348       cearray->AddLast(driftCalib);
2349     }
2350   }
2351   //
2352   // 0. Filter almost emty graphs
2353   //
2354
2355   for (Int_t i=0; i<72;i++){
2356     TGraph *graph= (TGraph*)arrT->At(i);
2357     if (!graph) continue;
2358     if (graph->GetN()<kMinPoints){
2359       arrT->AddAt(0,i);
2360       delete graph;  // delete empty graph
2361       continue;
2362     }
2363     if (tmin<0) tmin = graph->GetX()[0];
2364     if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2365     //
2366     if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2367     if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2368   }
2369   //
2370   // 1. calculate median and RMS per side
2371   //
2372   TArrayF arrA(100000), arrC(100000);
2373   Int_t nA=0, nC=0;
2374   Double_t medianA=0, medianC=0;
2375   Double_t rmsA=0, rmsC=0;
2376   for (Int_t isec=0; isec<72;isec++){
2377     TGraph *graph= (TGraph*)arrT->At(isec);
2378     if (!graph) continue;
2379     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2380       if (graph->GetY()[ipoint]<kMinTime) continue;
2381       if (nA>=arrA.fN) arrA.Set(nA*2);
2382       if (nC>=arrC.fN) arrC.Set(nC*2);
2383       if (isec%36<18)  arrA[nA++]= graph->GetY()[ipoint];
2384       if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2385     }
2386   }
2387   if (nA>0){
2388     medianA=TMath::Median(nA,arrA.fArray);
2389     rmsA   =TMath::RMS(nA,arrA.fArray);
2390   }
2391   if (nC>0){
2392     medianC=TMath::Median(nC,arrC.fArray);
2393     rmsC   =TMath::RMS(nC,arrC.fArray);
2394   }
2395   //
2396   // 2. Filter graphs - in respect with side medians
2397   //
2398   TArrayD vecX(100000), vecY(100000);
2399   for (Int_t isec=0; isec<72;isec++){
2400     TGraph *graph= (TGraph*)arrT->At(isec);
2401     if (!graph) continue;
2402     Double_t median = (isec%36<18) ? medianA: medianC;
2403     Double_t rms    = (isec%36<18) ? rmsA:    rmsC;
2404     Int_t naccept=0;
2405     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2406       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2407       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2408       vecX[naccept]= graph->GetX()[ipoint];
2409       vecY[naccept]= graph->GetY()[ipoint];
2410       naccept++;
2411     }
2412     if (naccept<kMinPoints){
2413       arrT->AddAt(0,isec);
2414       delete graph;  // delete empty graph
2415       continue;
2416     }
2417     TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2418     delete graph;
2419     arrT->AddAt(graph2,isec);
2420   }
2421   //
2422   // 3. Cut in respect wit the graph median
2423   //
2424   for (Int_t i=0; i<72;i++){
2425     TGraph *graph= (TGraph*)arrT->At(i);
2426     if (!graph) continue;
2427     //
2428     // filter in range
2429     //
2430     TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2431     if (!graphTS0) continue;
2432     if (graphTS0->GetN()<kMinPoints) {
2433       delete graphTS0;
2434       delete graph;
2435       arrT->AddAt(0,i);
2436       continue;
2437     }
2438     TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2439     graphTS->Sort();
2440     AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2441     if (pcstream){
2442       Int_t run = AliTPCcalibDB::Instance()->GetRun();
2443       (*pcstream)<<"filterCE"<<
2444         "run="<<run<<
2445         "isec="<<i<<
2446         "mY="<<medianY<<
2447         "graph.="<<graph<<
2448         "graphTS0.="<<graphTS0<<
2449         "graphTS.="<<graphTS<<
2450         "\n";
2451     }
2452     delete graphTS0;
2453     if (!graphTS) continue;
2454     arrT->AddAt(graphTS,i);
2455     delete graph;
2456   }
2457   //
2458   // Recalculate the mean time A side C side
2459   //
2460   TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2461   Int_t meanPoints=(nA+nC)/72;  // mean number of points
2462   for (Int_t itime=0; itime<200; itime++){
2463     nA=0, nC=0;
2464     Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2465     for (Int_t i=0; i<72;i++){
2466       TGraph *graph= (TGraph*)arrT->At(i);
2467       if (!graph) continue;
2468       if (graph->GetN()<(meanPoints/4)) continue;
2469       if ( (i%36)<18 )  arrA[nA++]=graph->Eval(time);
2470       if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2471     }
2472     xA[itime]=time;
2473     xC[itime]=time;
2474     yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2475     yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2476     eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2477     eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2478   }
2479   //
2480   Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2481   Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2482   if (pcstream){
2483     Int_t run = AliTPCcalibDB::Instance()->GetRun();
2484     (*pcstream)<<"filterAC"<<
2485       "run="<<run<<
2486       "nA="<<nA<<
2487       "nC="<<nC<<
2488       "rmsTA="<<rmsTA<<
2489       "rmsTC="<<rmsTC<<
2490       "\n";
2491   }
2492   //
2493   TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2494   TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2495   TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2496   if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2497   TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2498   if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2499   delete grA;
2500   delete grC;
2501   if (nA<kMinSectors) arrT->AddAt(0,72);
2502   else arrT->AddAt(graphTSA,72);
2503   if (nC<kMinSectors) arrT->AddAt(0,73);
2504   else arrT->AddAt(graphTSC,73);
2505 }
2506
2507
2508 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2509   //
2510   // Filter Drift velocity measurement using the tracks
2511   // 0.  remove outlyers - error based
2512   //     cutSigma
2513   //
2514   //
2515   const Int_t kMinPoints=1;  // minimal number of points to define value
2516   TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2517   Double_t medianY=0;
2518   if (!arrT) return;
2519   for (Int_t i=0; i<arrT->GetEntries();i++){
2520     TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2521     if (!graph) continue;
2522     if (graph->GetN()<kMinPoints){
2523       delete graph;
2524       arrT->AddAt(0,i);
2525       continue;
2526     }
2527     TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2528     if (!graph2) {
2529       delete graph; arrT->AddAt(0,i); continue;
2530     }
2531     if (graph2->GetN()<1) {
2532       delete graph; arrT->AddAt(0,i); continue;
2533     }
2534     graph2->SetName(graph->GetName());
2535     graph2->SetTitle(graph->GetTitle());
2536     arrT->AddAt(graph2,i);
2537     if (pcstream){
2538       (*pcstream)<<"filterTracks"<<
2539         "run="<<run<<
2540         "isec="<<i<<
2541         "mY="<<medianY<<
2542         "graph.="<<graph<<
2543         "graph2.="<<graph2<<
2544         "\n";
2545     }
2546     delete graph;
2547   }
2548 }
2549
2550
2551
2552
2553
2554 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2555   //
2556   //
2557   // get laser time offset
2558   // median around timeStamp+-deltaT
2559   // QA - chi2 needed for later usage - to be added
2560   //    - currently cut on error
2561   //
2562   Int_t kMinPoints=1;
2563   Double_t kMinDelay=0.01;
2564   Double_t kMinDelayErr=0.0001;
2565
2566   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2567   if (!array) return 0;
2568   TGraphErrors *tlaser=0;
2569   if (array){
2570     if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2571     if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2572   }
2573   if (!tlaser) return 0;
2574   Int_t npoints0= tlaser->GetN();
2575   if (npoints0==0) return 0;
2576   Double_t *xlaser = new Double_t[npoints0];
2577   Double_t *ylaser = new Double_t[npoints0];
2578   Int_t npoints=0;
2579   for (Int_t i=0;i<npoints0;i++){
2580     //printf("%d\n",i);
2581     if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2582     if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2583     xlaser[npoints]=tlaser->GetX()[npoints];
2584     ylaser[npoints]=tlaser->GetY()[npoints];
2585     npoints++;
2586   }
2587   //
2588   //
2589   Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2590   Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2591   //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2592   if (index0<0) index0=0;
2593   if (index1>=npoints-1) index1=npoints-1;
2594   if (index1-index0<kMinPoints) return 0;
2595   //
2596   //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2597     Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2598   delete [] ylaser;
2599   delete [] xlaser;
2600   return mean;
2601 }
2602
2603
2604
2605
2606 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2607   //
2608   // Filter Goofie data
2609   // goofieArray - points will be filtered
2610   // deltaT      - smmothing time window
2611   // cutSigma    - outler sigma cut in rms
2612   // minVn, maxVd- range absolute cut for variable vd/pt
2613   //             - to be tuned
2614   //
2615   // Ignore goofie if not enough points
2616   //
2617   const Int_t kMinPoints = 3;
2618   //
2619
2620   TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2621   TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2622   TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2623   TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2624   if (!graphvd) return;
2625   if (graphvd->GetN()<kMinPoints){
2626     delete graphvd;
2627     goofieArray->GetSensorNum(2)->SetGraph(0);
2628     return;
2629   }
2630   //
2631   // 1. Caluclate medians of critical variables
2632   //    drift velcocity
2633   //    P/T
2634   //    area near peak
2635   //    area far  peak
2636   //
2637   Double_t medianpt=0;
2638   Double_t medianvd=0, sigmavd=0;
2639   Double_t medianan=0;
2640   Double_t medianaf=0;
2641   Int_t    entries=graphvd->GetN();
2642   Double_t yvdn[10000];
2643   Int_t nvd=0;
2644   //
2645   for (Int_t ipoint=0; ipoint<entries; ipoint++){
2646     if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2647     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2648     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2649     yvdn[nvd++]=graphvd->GetY()[ipoint];
2650   }
2651   if (nvd<kMinPoints){
2652     delete graphvd;
2653     goofieArray->GetSensorNum(2)->SetGraph(0);
2654     return;
2655   }
2656   //
2657   Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2658   if (nuni>=kMinPoints){
2659     AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2660   }else{
2661     medianvd = TMath::Median(nvd, yvdn);
2662   }
2663
2664   TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2665   TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2666   TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2667   TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2668   TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2669   TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2670   delete graphpt0;
2671   delete graphpt1;
2672   delete graphan0;
2673   delete graphan1;
2674   delete graphaf0;
2675   delete graphaf1;
2676   //
2677   // 2. Make outlyer graph
2678   //
2679   Int_t nOK=0;
2680   TGraph graphOut(*graphvd);
2681   for (Int_t i=0; i<entries;i++){
2682     //
2683     Bool_t isOut=kFALSE;
2684     if (graphpt->GetY()[i]<=0.0000001) {  graphOut.GetY()[i]=1; continue;}
2685     if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) {  graphOut.GetY()[i]=1; continue;}
2686
2687     if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2688       isOut|=kTRUE;
2689     if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2690     if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2691     if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2692     graphOut.GetY()[i]= (isOut)?1:0;
2693     if (!isOut) nOK++;
2694   }
2695   if (nOK<kMinPoints) {
2696     delete graphvd;
2697     goofieArray->GetSensorNum(2)->SetGraph(0);
2698     return;
2699   }
2700   //
2701   // 3. Filter out outlyers - and smooth
2702   //
2703   TVectorF vmedianArray(goofieArray->NumSensors());
2704   TVectorF vrmsArray(goofieArray->NumSensors());
2705   Double_t xnew[10000];
2706   Double_t ynew[10000];
2707   TObjArray junk;
2708   junk.SetOwner(kTRUE);
2709   Bool_t isOK=kTRUE;
2710   //
2711   //
2712   for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2713     isOK=kTRUE;
2714     AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2715     TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2716     //
2717     if (!sensor) continue;
2718     graphOld = sensor->GetGraph();
2719     if (graphOld) {
2720       sensor->SetGraph(0);
2721       Int_t nused=0;
2722       for (Int_t i=0;i<entries;i++){
2723         if (graphOut.GetY()[i]>0.5) continue;
2724         xnew[nused]=graphOld->GetX()[i];
2725         ynew[nused]=graphOld->GetY()[i];
2726         nused++;
2727       }
2728       graphNew = new TGraph(nused,xnew,ynew);
2729       junk.AddLast(graphNew);
2730       junk.AddLast(graphOld);
2731       Double_t median=0;
2732       graphNew0  = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2733       if (graphNew0!=0){
2734         junk.AddLast(graphNew0);
2735         graphNew1  = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2736         if (graphNew1!=0){
2737           junk.AddLast(graphNew1);
2738           graphNew2  = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2739           if (graphNew2!=0) {
2740             vrmsArray[isensor]   =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2741             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2742             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2743             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2744             printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2745             vmedianArray[isensor]=median;
2746             //
2747           }
2748         }
2749       }
2750     }
2751     if (!graphOld) {  isOK=kFALSE; graphOld =&graphOut;}
2752     if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2753     if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2754     if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2755     (*pcstream)<<"goofieA"<<
2756       Form("isOK_%d.=",isensor)<<isOK<<
2757       Form("s_%d.=",isensor)<<sensor<<
2758       Form("gr_%d.=",isensor)<<graphOld<<
2759       Form("gr0_%d.=",isensor)<<graphNew0<<
2760       Form("gr1_%d.=",isensor)<<graphNew1<<
2761       Form("gr2_%d.=",isensor)<<graphNew2;
2762     if (isOK) sensor->SetGraph(graphNew2);
2763   }
2764   (*pcstream)<<"goofieA"<<
2765     "vmed.="<<&vmedianArray<<
2766     "vrms.="<<&vrmsArray<<
2767     "\n";
2768   junk.Delete();   // delete temoprary graphs
2769
2770 }
2771
2772
2773
2774
2775
2776 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2777   //
2778   // Make a statistic matrix
2779   // Input parameters:
2780   //   array        - TObjArray of AliRelKalmanAlign
2781   //   minFraction  - minimal ration of accepted tracks
2782   //   minStat      - minimal statistic (number of accepted tracks)
2783   //   maxvd        - maximal deviation for the 1
2784   // Output matrix:
2785   //    columns    - Mean, Median, RMS
2786   //    row        - parameter type (rotation[3], translation[3], drift[3])
2787   if (!array) return 0;
2788   if (array->GetEntries()<=0) return 0;
2789   //  Int_t entries = array->GetEntries();
2790   Int_t entriesFast = array->GetEntriesFast();
2791   TVectorD state(9);
2792   TVectorD *valArray[9];
2793   for (Int_t i=0; i<9; i++){
2794     valArray[i] = new TVectorD(entriesFast);
2795   }
2796   Int_t naccept=0;
2797   for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
2798     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
2799     if (!kalman) continue;
2800     if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
2801     if (kalman->GetNUpdates()<minStat) continue;
2802     if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
2803     kalman->GetState(state);
2804     for (Int_t ipar=0; ipar<9; ipar++)
2805       (*valArray[ipar])[naccept]=state[ipar];
2806     naccept++;
2807   }
2808   if (naccept<2) return 0;
2809   TMatrixD *pstat=new TMatrixD(9,3);
2810   TMatrixD &stat=*pstat;
2811   for (Int_t ipar=0; ipar<9; ipar++){
2812     stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
2813     stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
2814     stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
2815   }
2816   return pstat;
2817 }
2818
2819
2820 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
2821   //
2822   // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
2823   // Input:
2824   //   array     - input array
2825   //   stat      - mean parameters statistic
2826   //   direction -
2827   //   sigmaCut  - maximal allowed deviation from mean in terms of RMS
2828   if (!array) return 0;
2829   if (array->GetEntries()<=0) return 0;
2830   if (!(&stat)) return 0;
2831   // error increase in 1 hour
2832   const Double_t kerrsTime[9]={
2833     0.00001,  0.00001, 0.00001,
2834     0.001,    0.001,   0.001,
2835     0.0001,  0.001,   0.0001};
2836   //
2837   //
2838   Int_t entries = array->GetEntriesFast();
2839   TObjArray *sArray= new TObjArray(entries);
2840   AliRelAlignerKalman * sKalman =0;
2841   TVectorD state(9);
2842   for (Int_t i=0; i<entries; i++){
2843     Int_t index=(direction)? entries-i-1:i;
2844     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
2845     if (!kalman) continue;
2846     Bool_t isOK=kTRUE;
2847     kalman->GetState(state);
2848     for (Int_t ipar=0; ipar<9; ipar++){
2849       if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
2850     }
2851     if (!sKalman &&isOK) {
2852       sKalman=new AliRelAlignerKalman(*kalman);
2853       sKalman->SetRejectOutliers(kFALSE);
2854       sKalman->SetRunNumber(kalman->GetRunNumber());
2855       sKalman->SetTimeStamp(kalman->GetTimeStamp());
2856     }
2857     if (!sKalman) continue;
2858     Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
2859     for (Int_t ipar=0; ipar<9; ipar++){
2860 //       (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
2861 //       (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
2862 //       (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
2863       (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
2864     }
2865     sKalman->SetRunNumber(kalman->GetRunNumber());
2866     if (!isOK) sKalman->SetRunNumber(0);
2867     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2868     if (!isOK) continue;
2869     sKalman->SetRejectOutliers(kFALSE);
2870     sKalman->SetRunNumber(kalman->GetRunNumber());
2871     sKalman->SetTimeStamp(kalman->GetTimeStamp());
2872     sKalman->Merge(kalman);
2873     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2874     //sKalman->Print();
2875   }
2876   return sArray;
2877 }
2878
2879 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
2880   //
2881   // Merge 2 RelKalman arrays
2882   // Input:
2883   //   arrayP    - rel kalman in direction plus
2884   //   arrayM    - rel kalman in direction minus
2885   if (!arrayP) return 0;
2886   if (arrayP->GetEntries()<=0) return 0;
2887   if (!arrayM) return 0;
2888   if (arrayM->GetEntries()<=0) return 0;
2889   Int_t entries = arrayP->GetEntriesFast();
2890   TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
2891   for (Int_t i=0; i<entries; i++){
2892      AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
2893      AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
2894      if (!kalmanP) continue;
2895      if (!kalmanM) continue;
2896      AliRelAlignerKalman *kalman = new AliRelAlignerKalman(*kalmanP);
2897      kalman->Merge(kalmanM);
2898      array->AddAt(kalman,i);
2899   }
2900   return array;
2901 }