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