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 *const 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 Int_t AliTPCcalibDButil::GetCurrentReferenceRun(const char* type) const {
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 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   dx = 0;
1611   y = 0;
1612
1613   if(!graph) return 0;
1614   if(graph->GetN() < 1) return 0;
1615
1616   Int_t index=0;
1617   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1618   if (index<0) index=0;
1619   if(graph->GetN()==1) {
1620     dx = xref-graph->GetX()[index];
1621   }
1622   else {
1623     if (index>=graph->GetN()-1) index=graph->GetN()-2;
1624     if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1625     dx = xref-graph->GetX()[index];
1626   }
1627   y  = graph->GetY()[index];
1628   return index;
1629 }
1630
1631 Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1632   //
1633   // Get the correction of the trigger offset
1634   // combining information from the laser track calibration
1635   // and from cosmic calibration
1636   //
1637   // run       - run number
1638   // timeStamp - tim stamp in seconds
1639   // deltaT    - integration period to calculate offset
1640   // deltaTLaser -max validity of laser data
1641   // valType   - 0 - median, 1- mean
1642   //
1643   // Integration vaues are just recomendation - if not possible to get points
1644   // automatically increase the validity by factor 2
1645   // (recursive algorithm until one month of data taking)
1646   //
1647   //
1648   const Float_t kLaserCut=0.0005;
1649   const Int_t   kMaxPeriod=3600*24*30*12; // one year max
1650   const Int_t   kMinPoints=20;
1651   //
1652   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1653   if (!array) {
1654     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1655   }
1656   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1657   if (!array) return 0;
1658   //
1659   TGraphErrors *laserA[3]={0,0,0};
1660   TGraphErrors *laserC[3]={0,0,0};
1661   TGraphErrors *cosmicAll=0;
1662   laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1663   laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1664   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1665   //
1666   //
1667   if (!cosmicAll) return 0;
1668   Int_t nmeasC=cosmicAll->GetN();
1669   Float_t *tdelta = new Float_t[nmeasC];
1670   Int_t nused=0;
1671   for (Int_t i=0;i<nmeasC;i++){
1672     if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1673     Float_t ccosmic=cosmicAll->GetY()[i];
1674     Double_t yA=0,yC=0,dA=0,dC=0;
1675     if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1676     if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1677     //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1678     //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1679     //
1680     if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1681     Float_t claser=0;
1682     if (TMath::Abs(yA-yC)<kLaserCut) {
1683       claser=(yA-yC)*0.5;
1684     }else{
1685       if (i%2==0)  claser=yA;
1686       if (i%2==1)  claser=yC;
1687     }
1688     tdelta[nused]=ccosmic-claser;
1689     nused++;
1690   }
1691   if (nused<kMinPoints &&deltaT<kMaxPeriod) return  AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1692   if (nused<kMinPoints) {
1693     printf("AliFatal: No time offset calibration available\n");
1694     return 0;
1695   }
1696   Double_t median = TMath::Median(nused,tdelta);
1697   Double_t mean  = TMath::Mean(nused,tdelta);
1698   delete [] tdelta;
1699   return (valType==0) ? median:mean;
1700 }
1701
1702 Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1703   //
1704   // Get the correction of the drift velocity
1705   // combining information from the laser track calibration
1706   // and from cosmic calibration
1707   //
1708   // dist      - return value - distance to closest point in graph
1709   // run       - run number
1710   // timeStamp - tim stamp in seconds
1711   // deltaT    - integration period to calculate time0 offset
1712   // deltaTLaser -max validity of laser data
1713   // valType   - 0 - median, 1- mean
1714   //
1715   // Integration vaues are just recomendation - if not possible to get points
1716   // automatically increase the validity by factor 2
1717   // (recursive algorithm until one month of data taking)
1718   //
1719   //
1720   //
1721   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1722   if (!array) {
1723     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1724   }
1725   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1726   if (!array) return 0;
1727   TGraphErrors *cosmicAll=0;
1728   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1729   if (!cosmicAll) return 0;
1730   Double_t grY=0;
1731   AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1732
1733   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1734   Double_t vcosmic =  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1735   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1736   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
1737   return  vcosmic-t0;
1738
1739   /*
1740     Example usage:
1741
1742     Int_t run=89000
1743     TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1744     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1745     laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1746     //
1747     Double_t *yvd= new Double_t[cosmicAll->GetN()];
1748     Double_t *yt0= new Double_t[cosmicAll->GetN()];
1749     for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1750     for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1751
1752     TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1753     TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1754
1755   */
1756
1757 }
1758
1759 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1760 {
1761   //
1762   // Create a default name for the gui file
1763   //
1764
1765   return Form("guiRefTreeRun%s.root",GetRefValidity());
1766 }
1767
1768 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1769 {
1770   //
1771   // Create a gui reference tree
1772   // if dirname and filename are empty default values will be used
1773   // this is the recommended way of using this function
1774   // it allows to check whether a file with the given run validity alredy exists
1775   //
1776   if (!AliCDBManager::Instance()->GetDefaultStorage()){
1777     AliError("Default Storage not set. Cannot create reference calibration Tree!");
1778     return kFALSE;
1779   }
1780
1781   TString file=filename;
1782   if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1783
1784   AliTPCPreprocessorOnline prep;
1785   //noise and pedestals
1786   if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1787   if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1788   if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1789   //pulser data
1790   if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1791   if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1792   if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1793   if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1794   //CE data
1795   if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1796   if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1797   if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1798   if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1799   //Altro data
1800   if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1801   if (fRefALTROZsThr    ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr    )));
1802   if (fRefALTROFPED     ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED     )));
1803   if (fRefALTROAcqStop  ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop  )));
1804   if (fRefALTROMasked   ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked   )));
1805   //QA
1806   AliTPCdataQA *dataQA=fRefDataQA;
1807   if (dataQA) {
1808     if (dataQA->GetNLocalMaxima())
1809       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1810     if (dataQA->GetMaxCharge())
1811       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1812     if (dataQA->GetMeanCharge())
1813       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1814     if (dataQA->GetNoThreshold())
1815       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1816     if (dataQA->GetNTimeBins())
1817       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1818     if (dataQA->GetNPads())
1819       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1820     if (dataQA->GetTimePosition())
1821       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1822   }
1823   prep.DumpToFile(file.Data());
1824   return kTRUE;
1825 }
1826
1827 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1828   //
1829   // Get the correction of the drift velocity using the laser tracks calbration
1830   //
1831   // run       - run number
1832   // timeStamp - tim stamp in seconds
1833   // deltaT    - integration period to calculate time0 offset
1834   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1835   // Note in case no data form both A and C side - the value from active side used
1836   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1837   TGraphErrors *grlaserA=0;
1838   TGraphErrors *grlaserC=0;
1839   Double_t vlaserA=0, vlaserC=0;
1840   if (!array) return 0;
1841   grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1842   grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1843   Double_t deltaY;
1844   if (grlaserA && grlaserA->GetN()>0) {
1845     AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1846     if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
1847     else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1848   }
1849   if (grlaserC && grlaserC->GetN()>0) {
1850     AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1851     if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
1852     else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1853   }
1854   if (side==0) return vlaserA;
1855   if (side==1) return vlaserC;
1856   Double_t mdrift=(vlaserA+vlaserC)*0.5;
1857   if (!grlaserA) return vlaserC;
1858   if (!grlaserC) return vlaserA;
1859   return mdrift;
1860 }
1861
1862
1863 Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1864   //
1865   // Get the correction of the drift velocity using the CE laser data
1866   // combining information from the CE,  laser track calibration
1867   // and P/T calibration
1868   //
1869   // run       - run number
1870   // timeStamp - tim stamp in seconds
1871   // deltaT    - integration period to calculate time0 offset
1872   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1873   TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
1874   if (!arrT) return 0;
1875   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
1876   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1877   AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1878   //
1879   //
1880   Double_t corrPTA = 0, corrPTC=0;
1881   Double_t ltime0A = 0, ltime0C=0;
1882   Double_t gry=0;
1883   Double_t corrA=0, corrC=0;
1884   Double_t timeA=0, timeC=0;
1885   const Double_t kEpsilon = 0.00001;
1886   TGraph *graphA = (TGraph*)arrT->At(72);
1887   TGraph *graphC = (TGraph*)arrT->At(73);
1888   if (!graphA && !graphC) return 0.;
1889   if (graphA &&graphA->GetN()>0) {
1890     AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1891     timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1892     Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1893     ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1894     if(ltime0A < kEpsilon) return 0;
1895     if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
1896     corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1897     corrA-=corrPTA;
1898   }
1899   if (graphC&&graphC->GetN()>0){
1900     AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1901     timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1902     Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1903     ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1904     if(ltime0C < kEpsilon) return 0;
1905 if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
1906     corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1907     corrC-=corrPTC;
1908   }
1909
1910   if (side ==0 ) return corrA;
1911   if (side ==1 ) return corrC;
1912   Double_t corrM= (corrA+corrC)*0.5;
1913   if (!graphA) corrM=corrC;
1914   if (!graphC) corrM=corrA;
1915   return corrM;
1916 }
1917
1918 Double_t  AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1919   //
1920   // return drift velocity using the TPC-ITS matchin method
1921   // return also distance to the closest point
1922   //
1923   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1924   TGraphErrors *graph=0;
1925   dist=0;
1926   if (!array) return 0;
1927   //array->ls();
1928   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
1929   if (!graph) return 0;
1930   Double_t deltaY;
1931   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1932   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1933   return value;
1934 }
1935
1936 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
1937   //
1938   // Get time dependent time 0 (trigger delay in cm) correction
1939   // Arguments:
1940   // timestamp - timestamp
1941   // run       - run number
1942   //
1943   // Notice - Extrapolation outside of calibration range  - using constant function
1944   //
1945   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1946   TGraphErrors *graph=0;
1947   dist=0;
1948   if (!array) return 0;
1949   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
1950   if (!graph) return 0;
1951   Double_t deltaY;
1952   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
1953   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
1954   return value;
1955 }
1956
1957
1958
1959
1960
1961 Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1962   //
1963   // VERY obscure method - we need something in framework
1964   // Find the TPC runs with temperature OCDB entry
1965   // cache the start and end of the run
1966   //
1967   AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1968   if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1969   if (!storage) return 0;
1970   TString path=storage->GetURI();
1971   TString runsT;
1972   {
1973     TString command;
1974     if (path.Contains("local")){  // find the list if local system
1975       path.ReplaceAll("local://","");
1976       path+="TPC/Calib/Temperature";
1977       command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
1978     }
1979     runsT=gSystem->GetFromPipe(command);
1980   }
1981   TObjArray *arr= runsT.Tokenize("r");
1982   if (!arr) return 0;
1983   //
1984   TArrayI indexes(arr->GetEntries());
1985   TArrayI runs(arr->GetEntries());
1986   Int_t naccept=0;
1987   {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1988       Int_t irunN = atoi(arr->At(irun)->GetName());
1989       if (irunN<startRun) continue;
1990       if (irunN>stopRun) continue;
1991       runs[naccept]=irunN;
1992       naccept++;
1993     }}
1994   fRuns.Set(naccept);
1995   fRunsStart.Set(fRuns.fN);
1996   fRunsStop.Set(fRuns.fN);
1997   TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1998   for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
1999
2000   //
2001   AliCDBEntry * entry = 0;
2002   {for (Int_t irun=0;irun<fRuns.fN; irun++){
2003       entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2004       if (!entry) continue;
2005       AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2006       if (!tmpRun) continue;
2007       fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2008       fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2009       //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2010     }}
2011   return fRuns.fN;
2012 }
2013
2014
2015 Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
2016   //
2017   // binary search - find the run for given time stamp
2018   //
2019   Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2020   Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2021   Int_t cindex  = -1;
2022   for (Int_t index=index0; index<=index1; index++){
2023     if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2024     if (debug) {
2025       printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2026     }
2027   }
2028   if (cindex<0) cindex =(index0+index1)/2;
2029   if (cindex<0) {
2030     return 0;
2031   }
2032   return fRuns[cindex];
2033 }
2034
2035
2036
2037
2038
2039 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2040   //
2041   // filter outlyer measurement
2042   // Only points around median +- sigmaCut filtered
2043   if (!graph) return  0;
2044   Int_t kMinPoints=2;
2045   Int_t npoints0 = graph->GetN();
2046   Int_t npoints=0;
2047   Float_t  rmsY=0;
2048   //
2049   //
2050   if (npoints0<kMinPoints) return 0;
2051
2052   Double_t *outx=new Double_t[npoints0];
2053   Double_t *outy=new Double_t[npoints0];
2054   for (Int_t iter=0; iter<3; iter++){
2055     npoints=0;
2056     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2057       if (graph->GetY()[ipoint]==0) continue;
2058       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2059       outx[npoints]  = graph->GetX()[ipoint];
2060       outy[npoints]  = graph->GetY()[ipoint];
2061       npoints++;
2062     }
2063     if (npoints<=1) break;
2064     medianY  =TMath::Median(npoints,outy);
2065     rmsY   =TMath::RMS(npoints,outy);
2066   }
2067   TGraph *graphOut=0;
2068   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2069   delete [] outx;
2070   delete [] outy;
2071   return graphOut;
2072 }
2073
2074
2075 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2076   //
2077   // filter outlyer measurement
2078   // Only points around median +- cut filtered
2079   if (!graph) return  0;
2080   Int_t kMinPoints=2;
2081   Int_t npoints0 = graph->GetN();
2082   Int_t npoints=0;
2083   Float_t  rmsY=0;
2084   //
2085   //
2086   if (npoints0<kMinPoints) return 0;
2087
2088   Double_t *outx=new Double_t[npoints0];
2089   Double_t *outy=new Double_t[npoints0];
2090   for (Int_t iter=0; iter<3; iter++){
2091     npoints=0;
2092     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2093       if (graph->GetY()[ipoint]==0) continue;
2094       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2095       outx[npoints]  = graph->GetX()[ipoint];
2096       outy[npoints]  = graph->GetY()[ipoint];
2097       npoints++;
2098     }
2099     if (npoints<=1) break;
2100     medianY  =TMath::Median(npoints,outy);
2101     rmsY   =TMath::RMS(npoints,outy);
2102   }
2103   TGraph *graphOut=0;
2104   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2105   delete [] outx;
2106   delete [] outy;
2107   return graphOut;
2108 }
2109
2110
2111
2112 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2113   //
2114   // filter outlyer measurement
2115   // Only points with normalized errors median +- sigmaCut filtered
2116   //
2117   Int_t kMinPoints=10;
2118   Int_t npoints0 = graph->GetN();
2119   Int_t npoints=0;
2120   Float_t  medianErr=0, rmsErr=0;
2121   //
2122   //
2123   if (npoints0<kMinPoints) return 0;
2124
2125   Double_t *outx=new Double_t[npoints0];
2126   Double_t *outy=new Double_t[npoints0];
2127   Double_t *erry=new Double_t[npoints0];
2128   Double_t *nerry=new Double_t[npoints0];
2129   Double_t *errx=new Double_t[npoints0];
2130
2131   for (Int_t iter=0; iter<3; iter++){
2132     npoints=0;
2133     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2134       nerry[npoints]  = graph->GetErrorY(ipoint);
2135       if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2136       erry[npoints]  = graph->GetErrorY(ipoint);
2137       outx[npoints]  = graph->GetX()[ipoint];
2138       outy[npoints]  = graph->GetY()[ipoint];
2139       errx[npoints]  = graph->GetErrorY(ipoint);
2140       npoints++;
2141     }
2142     if (npoints==0) break;
2143     medianErr=TMath::Median(npoints,erry);
2144     medianY  =TMath::Median(npoints,outy);
2145     rmsErr   =TMath::RMS(npoints,erry);
2146   }
2147   TGraphErrors *graphOut=0;
2148   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2149   delete []outx;
2150   delete []outy;
2151   delete []erry;
2152   delete []nerry;
2153   delete []errx;
2154   return graphOut;
2155 }
2156
2157
2158 void AliTPCcalibDButil::Sort(TGraph *graph){
2159   //
2160   // sort array - neccessay for approx
2161   //
2162   Int_t npoints = graph->GetN();
2163   Int_t *indexes=new Int_t[npoints];
2164   Double_t *outx=new Double_t[npoints];
2165   Double_t *outy=new Double_t[npoints];
2166   TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2167   for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2168   for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2169   for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2170   for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2171
2172   delete [] indexes;
2173   delete [] outx;
2174   delete [] outy;
2175 }
2176 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2177   //
2178   // smmoth graph - mean on the interval
2179   //
2180   Sort(graph);
2181   Int_t npoints = graph->GetN();
2182   Double_t *outy=new Double_t[npoints];
2183
2184   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2185     Double_t lx=graph->GetX()[ipoint];
2186     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2187     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2188     if (index0<0) index0=0;
2189     if (index1>=npoints-1) index1=npoints-1;
2190     if ((index1-index0)>1){
2191       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2192     }else{
2193       outy[ipoint]=graph->GetY()[ipoint];
2194     }
2195   }
2196  //  TLinearFitter  fitter(3,"pol2");
2197 //   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2198 //     Double_t lx=graph->GetX()[ipoint];
2199 //     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2200 //     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2201 //     if (index0<0) index0=0;
2202 //     if (index1>=npoints-1) index1=npoints-1;
2203 //     fitter.ClearPoints();
2204 //     for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2205 //     if ((index1-index0)>1){
2206 //       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2207 //     }else{
2208 //       outy[ipoint]=graph->GetY()[ipoint];
2209 //     }
2210 //   }
2211
2212
2213
2214   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2215     graph->GetY()[ipoint] = outy[ipoint];
2216   }
2217   delete[] outy;
2218 }
2219
2220 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2221   //
2222   // Use constant interpolation outside of range
2223   //
2224   if (!graph) {
2225     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2226     return 0;
2227   }
2228
2229   if (graph->GetN()<1){
2230     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
2231     return 0;
2232   }
2233
2234
2235   if (xref<graph->GetX()[0]) return graph->GetY()[0];
2236   if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2237
2238   printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
2239
2240   if(graph->GetN()==1)
2241     return graph->Eval(graph->GetX()[0]);
2242
2243
2244   return graph->Eval(xref);
2245 }
2246
2247 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy,  Double_t sigmaCut){
2248   //
2249   // Filter DCS sensor information
2250   //   ymin     - minimal value
2251   //   ymax     - max value
2252   //   maxdy    - maximal deirivative
2253   //   sigmaCut - cut on values and derivative in terms of RMS distribution
2254   // Return value - accepted fraction
2255   //
2256   // Algorithm:
2257   //
2258   // 0. Calculate median and rms of values in specified range
2259   // 1. Filter out outliers - median+-sigmaCut*rms
2260   //    values replaced by median
2261   //
2262   AliSplineFit * fit    = sensor->GetFit();
2263   if (!fit) return 0.;
2264   Int_t          nknots = fit->GetKnots();
2265   if (nknots==0) {
2266     delete fit;
2267     sensor->SetFit(0);
2268     return 0;
2269   }
2270   //
2271   Double_t *yin0  = new Double_t[nknots];
2272   Double_t *yin1  = new Double_t[nknots];
2273   Int_t naccept=0;
2274
2275   for (Int_t iknot=0; iknot< nknots; iknot++){
2276     if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2277       yin0[naccept]  = fit->GetY0()[iknot];
2278       yin1[naccept]  = fit->GetY1()[iknot];
2279       if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2280       naccept++;
2281     }
2282   }
2283   if (naccept<1) {
2284     delete fit;
2285     sensor->SetFit(0);
2286     delete [] yin0;
2287     delete [] yin1;
2288     return 0.;
2289   }
2290
2291   Double_t medianY0=0, medianY1=0;
2292   Double_t rmsY0   =0, rmsY1=0;
2293   medianY0 = TMath::Median(naccept, yin0);
2294   medianY1 = TMath::Median(naccept, yin1);
2295   rmsY0    = TMath::RMS(naccept, yin0);
2296   rmsY1    = TMath::RMS(naccept, yin1);
2297   naccept=0;
2298   //
2299   // 1. Filter out outliers - median+-sigmaCut*rms
2300   //    values replaced by median
2301   //    if replaced the derivative set to 0
2302   //
2303   for (Int_t iknot=0; iknot< nknots; iknot++){
2304     Bool_t isOK=kTRUE;
2305     if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2306     if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2307     if (nknots<2) fit->GetY1()[iknot]=0;
2308     if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2309     if (!isOK){
2310       fit->GetY0()[iknot]=medianY0;
2311       fit->GetY1()[iknot]=0;
2312     }else{
2313       naccept++;
2314     }
2315   }
2316   delete [] yin0;
2317   delete [] yin1;
2318   return Float_t(naccept)/Float_t(nknots);
2319 }
2320
2321 Float_t  AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2322   //
2323   // Filter temperature array
2324   // tempArray    - array of temperatures         -
2325   // ymin         - minimal accepted temperature  - default 15
2326   // ymax         - maximal accepted temperature  - default 22
2327   // sigmaCut     - values filtered on interval median+-sigmaCut*rms - defaut 5
2328   // return value - fraction of filtered sensors
2329   const Double_t kMaxDy=0.1;
2330   Int_t nsensors=tempArray->NumSensors();
2331   if (nsensors==0) return 0.;
2332   Int_t naccept=0;
2333   for (Int_t isensor=0; isensor<nsensors; isensor++){
2334     AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2335     if (!sensor) continue;
2336     //printf("%d\n",isensor);
2337     FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2338     if (sensor->GetFit()==0){
2339       //delete sensor;
2340       tempArray->RemoveSensorNum(isensor);
2341     }else{
2342       naccept++;
2343     }
2344   }
2345   return Float_t(naccept)/Float_t(nsensors);
2346 }
2347
2348
2349 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2350   //
2351   // Filter CE data
2352   // Input parameters:
2353   //    deltaT   - smoothing window (in seconds)
2354   //    cutAbs   - max distance of the time info to the median (in time bins)
2355   //    cutSigma - max distance (in the RMS)
2356   //    pcstream - optional debug streamer to store original and filtered info
2357   // Hardwired parameters:
2358   //    kMinPoints =10;       // minimal number of points to define the CE
2359   //    kMinSectors=12;       // minimal number of sectors to define sideCE
2360   // Algorithm:
2361   // 0. Filter almost emty graphs (kMinPoints=10)
2362   // 1. calculate median and RMS per side
2363   // 2. Filter graphs - in respect with side medians
2364   //                  - cutAbs and cutDelta used
2365   // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2366   // 4. Calculate mean for A side and C side
2367   //
2368   const Int_t kMinPoints =10;       // minimal number of points to define the CE
2369   const Int_t kMinSectors=12;       // minimal number of sectors to define sideCE
2370   const Int_t kMinTime   =400;     // minimal arrival time of CE
2371   TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2372   Double_t medianY=0;
2373   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
2374   if (!cearray) return;
2375   Double_t tmin=-1;
2376   Double_t tmax=-1;
2377   //
2378   //
2379   AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2380   AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2381   if ( tempMapCE && cavernPressureCE){
2382     //
2383     //     Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2384     //     FilterSensor(cavernPressureCE,960,1050,10, 5.);
2385     //     if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2386     Bool_t isOK=kTRUE;
2387     if (isOK)  {
2388       // recalculate P/T correction map for time of the CE
2389       AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2390       driftCalib->SetName("driftPTCE");
2391       driftCalib->SetTitle("driftPTCE");
2392       cearray->AddLast(driftCalib);
2393     }
2394   }
2395   //
2396   // 0. Filter almost emty graphs
2397   //
2398
2399   for (Int_t i=0; i<72;i++){
2400     TGraph *graph= (TGraph*)arrT->At(i);
2401     if (!graph) continue;
2402     if (graph->GetN()<kMinPoints){
2403       arrT->AddAt(0,i);
2404       delete graph;  // delete empty graph
2405       continue;
2406     }
2407     if (tmin<0) tmin = graph->GetX()[0];
2408     if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2409     //
2410     if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2411     if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2412   }
2413   //
2414   // 1. calculate median and RMS per side
2415   //
2416   TArrayF arrA(100000), arrC(100000);
2417   Int_t nA=0, nC=0;
2418   Double_t medianA=0, medianC=0;
2419   Double_t rmsA=0, rmsC=0;
2420   for (Int_t isec=0; isec<72;isec++){
2421     TGraph *graph= (TGraph*)arrT->At(isec);
2422     if (!graph) continue;
2423     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2424       if (graph->GetY()[ipoint]<kMinTime) continue;
2425       if (nA>=arrA.fN) arrA.Set(nA*2);
2426       if (nC>=arrC.fN) arrC.Set(nC*2);
2427       if (isec%36<18)  arrA[nA++]= graph->GetY()[ipoint];
2428       if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2429     }
2430   }
2431   if (nA>0){
2432     medianA=TMath::Median(nA,arrA.fArray);
2433     rmsA   =TMath::RMS(nA,arrA.fArray);
2434   }
2435   if (nC>0){
2436     medianC=TMath::Median(nC,arrC.fArray);
2437     rmsC   =TMath::RMS(nC,arrC.fArray);
2438   }
2439   //
2440   // 2. Filter graphs - in respect with side medians
2441   //
2442   TArrayD vecX(100000), vecY(100000);
2443   for (Int_t isec=0; isec<72;isec++){
2444     TGraph *graph= (TGraph*)arrT->At(isec);
2445     if (!graph) continue;
2446     Double_t median = (isec%36<18) ? medianA: medianC;
2447     Double_t rms    = (isec%36<18) ? rmsA:    rmsC;
2448     Int_t naccept=0;
2449     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2450       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2451       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2452       vecX[naccept]= graph->GetX()[ipoint];
2453       vecY[naccept]= graph->GetY()[ipoint];
2454       naccept++;
2455     }
2456     if (naccept<kMinPoints){
2457       arrT->AddAt(0,isec);
2458       delete graph;  // delete empty graph
2459       continue;
2460     }
2461     TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2462     delete graph;
2463     arrT->AddAt(graph2,isec);
2464   }
2465   //
2466   // 3. Cut in respect wit the graph median
2467   //
2468   for (Int_t i=0; i<72;i++){
2469     TGraph *graph= (TGraph*)arrT->At(i);
2470     if (!graph) continue;
2471     //
2472     // filter in range
2473     //
2474     TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2475     if (!graphTS0) continue;
2476     if (graphTS0->GetN()<kMinPoints) {
2477       delete graphTS0;
2478       delete graph;
2479       arrT->AddAt(0,i);
2480       continue;
2481     }
2482     TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2483     graphTS->Sort();
2484     AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2485     if (pcstream){
2486       Int_t run = AliTPCcalibDB::Instance()->GetRun();
2487       (*pcstream)<<"filterCE"<<
2488         "run="<<run<<
2489         "isec="<<i<<
2490         "mY="<<medianY<<
2491         "graph.="<<graph<<
2492         "graphTS0.="<<graphTS0<<
2493         "graphTS.="<<graphTS<<
2494         "\n";
2495     }
2496     delete graphTS0;
2497     if (!graphTS) continue;
2498     arrT->AddAt(graphTS,i);
2499     delete graph;
2500   }
2501   //
2502   // Recalculate the mean time A side C side
2503   //
2504   TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2505   Int_t meanPoints=(nA+nC)/72;  // mean number of points
2506   for (Int_t itime=0; itime<200; itime++){
2507     nA=0, nC=0;
2508     Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2509     for (Int_t i=0; i<72;i++){
2510       TGraph *graph= (TGraph*)arrT->At(i);
2511       if (!graph) continue;
2512       if (graph->GetN()<(meanPoints/4)) continue;
2513       if ( (i%36)<18 )  arrA[nA++]=graph->Eval(time);
2514       if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2515     }
2516     xA[itime]=time;
2517     xC[itime]=time;
2518     yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2519     yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2520     eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2521     eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2522   }
2523   //
2524   Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2525   Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2526   if (pcstream){
2527     Int_t run = AliTPCcalibDB::Instance()->GetRun();
2528     (*pcstream)<<"filterAC"<<
2529       "run="<<run<<
2530       "nA="<<nA<<
2531       "nC="<<nC<<
2532       "rmsTA="<<rmsTA<<
2533       "rmsTC="<<rmsTC<<
2534       "\n";
2535   }
2536   //
2537   TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2538   TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2539   TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2540   if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2541   TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2542   if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2543   delete grA;
2544   delete grC;
2545   if (nA<kMinSectors) arrT->AddAt(0,72);
2546   else arrT->AddAt(graphTSA,72);
2547   if (nC<kMinSectors) arrT->AddAt(0,73);
2548   else arrT->AddAt(graphTSC,73);
2549 }
2550
2551
2552 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2553   //
2554   // Filter Drift velocity measurement using the tracks
2555   // 0.  remove outlyers - error based
2556   //     cutSigma
2557   //
2558   //
2559   const Int_t kMinPoints=1;  // minimal number of points to define value
2560   TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2561   Double_t medianY=0;
2562   if (!arrT) return;
2563   for (Int_t i=0; i<arrT->GetEntries();i++){
2564     TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
2565     if (!graph) continue;
2566     if (graph->GetN()<kMinPoints){
2567       delete graph;
2568       arrT->AddAt(0,i);
2569       continue;
2570     }
2571     TGraphErrors *graph2 = NULL;
2572     if(graph->GetN()<10) {
2573       graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2574       if (!graph2) {
2575         delete graph; arrT->AddAt(0,i); continue;
2576       }
2577     }
2578     else {
2579       graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2580       if (!graph2) {
2581         delete graph; arrT->AddAt(0,i); continue;
2582       }
2583     }
2584     if (graph2->GetN()<1) {
2585       delete graph; arrT->AddAt(0,i); continue;
2586     }
2587     graph2->SetName(graph->GetName());
2588     graph2->SetTitle(graph->GetTitle());
2589     arrT->AddAt(graph2,i);
2590     if (pcstream){
2591       (*pcstream)<<"filterTracks"<<
2592         "run="<<run<<
2593         "isec="<<i<<
2594         "mY="<<medianY<<
2595         "graph.="<<graph<<
2596         "graph2.="<<graph2<<
2597         "\n";
2598     }
2599     delete graph;
2600   }
2601 }
2602
2603
2604
2605
2606
2607 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2608   //
2609   //
2610   // get laser time offset
2611   // median around timeStamp+-deltaT
2612   // QA - chi2 needed for later usage - to be added
2613   //    - currently cut on error
2614   //
2615   Int_t kMinPoints=1;
2616   Double_t kMinDelay=0.01;
2617   Double_t kMinDelayErr=0.0001;
2618
2619   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2620   if (!array) return 0;
2621   TGraphErrors *tlaser=0;
2622   if (array){
2623     if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2624     if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2625   }
2626   if (!tlaser) return 0;
2627   Int_t npoints0= tlaser->GetN();
2628   if (npoints0==0) return 0;
2629   Double_t *xlaser = new Double_t[npoints0];
2630   Double_t *ylaser = new Double_t[npoints0];
2631   Int_t npoints=0;
2632   for (Int_t i=0;i<npoints0;i++){
2633     //printf("%d\n",i);
2634     if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2635     if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2636     xlaser[npoints]=tlaser->GetX()[npoints];
2637     ylaser[npoints]=tlaser->GetY()[npoints];
2638     npoints++;
2639   }
2640   //
2641   //
2642   Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2643   Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2644   //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2645   if (index0<0) index0=0;
2646   if (index1>=npoints-1) index1=npoints-1;
2647   if (index1-index0<kMinPoints) return 0;
2648   //
2649   //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2650     Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2651   delete [] ylaser;
2652   delete [] xlaser;
2653   return mean;
2654 }
2655
2656
2657
2658
2659 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2660   //
2661   // Filter Goofie data
2662   // goofieArray - points will be filtered
2663   // deltaT      - smmothing time window
2664   // cutSigma    - outler sigma cut in rms
2665   // minVn, maxVd- range absolute cut for variable vd/pt
2666   //             - to be tuned
2667   //
2668   // Ignore goofie if not enough points
2669   //
2670   const Int_t kMinPoints = 3;
2671   //
2672
2673   TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2674   TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2675   TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2676   TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2677   if (!graphvd) return;
2678   if (graphvd->GetN()<kMinPoints){
2679     delete graphvd;
2680     goofieArray->GetSensorNum(2)->SetGraph(0);
2681     return;
2682   }
2683   //
2684   // 1. Caluclate medians of critical variables
2685   //    drift velcocity
2686   //    P/T
2687   //    area near peak
2688   //    area far  peak
2689   //
2690   Double_t medianpt=0;
2691   Double_t medianvd=0, sigmavd=0;
2692   Double_t medianan=0;
2693   Double_t medianaf=0;
2694   Int_t    entries=graphvd->GetN();
2695   Double_t yvdn[10000];
2696   Int_t nvd=0;
2697   //
2698   for (Int_t ipoint=0; ipoint<entries; ipoint++){
2699     if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2700     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2701     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2702     yvdn[nvd++]=graphvd->GetY()[ipoint];
2703   }
2704   if (nvd<kMinPoints){
2705     delete graphvd;
2706     goofieArray->GetSensorNum(2)->SetGraph(0);
2707     return;
2708   }
2709   //
2710   Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2711   if (nuni>=kMinPoints){
2712     AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2713   }else{
2714     medianvd = TMath::Median(nvd, yvdn);
2715   }
2716
2717   TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2718   TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2719   TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2720   TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2721   TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2722   TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2723   delete graphpt0;
2724   delete graphpt1;
2725   delete graphan0;
2726   delete graphan1;
2727   delete graphaf0;
2728   delete graphaf1;
2729   //
2730   // 2. Make outlyer graph
2731   //
2732   Int_t nOK=0;
2733   TGraph graphOut(*graphvd);
2734   for (Int_t i=0; i<entries;i++){
2735     //
2736     Bool_t isOut=kFALSE;
2737     if (graphpt->GetY()[i]<=0.0000001) {  graphOut.GetY()[i]=1; continue;}
2738     if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) {  graphOut.GetY()[i]=1; continue;}
2739
2740     if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2741       isOut|=kTRUE;
2742     if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2743     if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2744     if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2745     graphOut.GetY()[i]= (isOut)?1:0;
2746     if (!isOut) nOK++;
2747   }
2748   if (nOK<kMinPoints) {
2749     delete graphvd;
2750     goofieArray->GetSensorNum(2)->SetGraph(0);
2751     return;
2752   }
2753   //
2754   // 3. Filter out outlyers - and smooth
2755   //
2756   TVectorF vmedianArray(goofieArray->NumSensors());
2757   TVectorF vrmsArray(goofieArray->NumSensors());
2758   Double_t xnew[10000];
2759   Double_t ynew[10000];
2760   TObjArray junk;
2761   junk.SetOwner(kTRUE);
2762   Bool_t isOK=kTRUE;
2763   //
2764   //
2765   for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2766     isOK=kTRUE;
2767     AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2768     TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2769     //
2770     if (!sensor) continue;
2771     graphOld = sensor->GetGraph();
2772     if (graphOld) {
2773       sensor->SetGraph(0);
2774       Int_t nused=0;
2775       for (Int_t i=0;i<entries;i++){
2776         if (graphOut.GetY()[i]>0.5) continue;
2777         xnew[nused]=graphOld->GetX()[i];
2778         ynew[nused]=graphOld->GetY()[i];
2779         nused++;
2780       }
2781       graphNew = new TGraph(nused,xnew,ynew);
2782       junk.AddLast(graphNew);
2783       junk.AddLast(graphOld);
2784       Double_t median=0;
2785       graphNew0  = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2786       if (graphNew0!=0){
2787         junk.AddLast(graphNew0);
2788         graphNew1  = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2789         if (graphNew1!=0){
2790           junk.AddLast(graphNew1);
2791           graphNew2  = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2792           if (graphNew2!=0) {
2793             vrmsArray[isensor]   =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2794             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2795             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2796             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2797             printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2798             vmedianArray[isensor]=median;
2799             //
2800           }
2801         }
2802       }
2803     }
2804     if (!graphOld) {  isOK=kFALSE; graphOld =&graphOut;}
2805     if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2806     if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2807     if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2808     (*pcstream)<<"goofieA"<<
2809       Form("isOK_%d.=",isensor)<<isOK<<
2810       Form("s_%d.=",isensor)<<sensor<<
2811       Form("gr_%d.=",isensor)<<graphOld<<
2812       Form("gr0_%d.=",isensor)<<graphNew0<<
2813       Form("gr1_%d.=",isensor)<<graphNew1<<
2814       Form("gr2_%d.=",isensor)<<graphNew2;
2815     if (isOK) sensor->SetGraph(graphNew2);
2816   }
2817   (*pcstream)<<"goofieA"<<
2818     "vmed.="<<&vmedianArray<<
2819     "vrms.="<<&vrmsArray<<
2820     "\n";
2821   junk.Delete();   // delete temoprary graphs
2822
2823 }
2824
2825
2826
2827
2828
2829 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2830   //
2831   // Make a statistic matrix
2832   // Input parameters:
2833   //   array        - TObjArray of AliRelKalmanAlign
2834   //   minFraction  - minimal ration of accepted tracks
2835   //   minStat      - minimal statistic (number of accepted tracks)
2836   //   maxvd        - maximal deviation for the 1
2837   // Output matrix:
2838   //    columns    - Mean, Median, RMS
2839   //    row        - parameter type (rotation[3], translation[3], drift[3])
2840   if (!array) return 0;
2841   if (array->GetEntries()<=0) return 0;
2842   //  Int_t entries = array->GetEntries();
2843   Int_t entriesFast = array->GetEntriesFast();
2844   TVectorD state(9);
2845   TVectorD *valArray[9];
2846   for (Int_t i=0; i<9; i++){
2847     valArray[i] = new TVectorD(entriesFast);
2848   }
2849   Int_t naccept=0;
2850   for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
2851     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
2852     if (!kalman) continue;
2853     if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
2854     if (kalman->GetNUpdates()<minStat) continue;
2855     if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
2856     kalman->GetState(state);
2857     for (Int_t ipar=0; ipar<9; ipar++)
2858       (*valArray[ipar])[naccept]=state[ipar];
2859     naccept++;
2860   }
2861   //if (naccept<2) return 0;
2862   if (naccept<1) return 0;
2863   TMatrixD *pstat=new TMatrixD(9,3);
2864   TMatrixD &stat=*pstat;
2865   for (Int_t ipar=0; ipar<9; ipar++){
2866     stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
2867     stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
2868     stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
2869   }
2870   return pstat;
2871 }
2872
2873
2874 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
2875   //
2876   // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
2877   // Input:
2878   //   array     - input array
2879   //   stat      - mean parameters statistic
2880   //   direction -
2881   //   sigmaCut  - maximal allowed deviation from mean in terms of RMS
2882   if (!array) return 0;
2883   if (array->GetEntries()<=0) return 0;
2884   if (!(&stat)) return 0;
2885   // error increase in 1 hour
2886   const Double_t kerrsTime[9]={
2887     0.00001,  0.00001, 0.00001,
2888     0.001,    0.001,   0.001,
2889     0.002,  0.01,   0.001};
2890   //
2891   //
2892   Int_t entries = array->GetEntriesFast();
2893   TObjArray *sArray= new TObjArray(entries);
2894   AliRelAlignerKalman * sKalman =0;
2895   TVectorD state(9);
2896   for (Int_t i=0; i<entries; i++){
2897     Int_t index=(direction)? entries-i-1:i;
2898     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
2899     if (!kalman) continue;
2900     Bool_t isOK=kTRUE;
2901     kalman->GetState(state);
2902     for (Int_t ipar=0; ipar<9; ipar++){
2903       if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
2904     }
2905     if (!sKalman &&isOK) {
2906       sKalman=new AliRelAlignerKalman(*kalman);
2907       sKalman->SetRejectOutliers(kFALSE);
2908       sKalman->SetRunNumber(kalman->GetRunNumber());
2909       sKalman->SetTimeStamp(kalman->GetTimeStamp());
2910     }
2911     if (!sKalman) continue;
2912     Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
2913     for (Int_t ipar=0; ipar<9; ipar++){
2914 //       (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
2915 //       (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
2916 //       (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
2917       (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
2918     }
2919     sKalman->SetRunNumber(kalman->GetRunNumber());
2920     if (!isOK) sKalman->SetRunNumber(0);
2921     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2922     if (!isOK) continue;
2923     sKalman->SetRejectOutliers(kFALSE);
2924     sKalman->SetRunNumber(kalman->GetRunNumber());
2925     sKalman->SetTimeStamp(kalman->GetTimeStamp());
2926     sKalman->Merge(kalman);
2927     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
2928     //sKalman->Print();
2929   }
2930   return sArray;
2931 }
2932
2933 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
2934   //
2935   // Merge 2 RelKalman arrays
2936   // Input:
2937   //   arrayP    - rel kalman in direction plus
2938   //   arrayM    - rel kalman in direction minus
2939   if (!arrayP) return 0;
2940   if (arrayP->GetEntries()<=0) return 0;
2941   if (!arrayM) return 0;
2942   if (arrayM->GetEntries()<=0) return 0;
2943   Int_t entries = arrayP->GetEntriesFast();
2944   TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
2945   for (Int_t i=0; i<entries; i++){
2946      AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
2947      AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
2948      if (!kalmanP) continue;
2949      if (!kalmanM) continue;
2950      AliRelAlignerKalman *kalman = new AliRelAlignerKalman(*kalmanP);
2951      kalman->Merge(kalmanM);
2952      array->AddAt(kalman,i);
2953   }
2954   return array;
2955 }