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