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 <AliCDBStorage.h>
  34 #include <AliDCSSensorArray.h>
  35 #include <AliTPCSensorTempArray.h>
  36 #include <AliDCSSensor.h>
  37 #include "AliTPCcalibDB.h"
  38 #include "AliTPCCalPad.h"
  39 #include "AliTPCCalROC.h"
  40 #include "AliTPCROC.h"
  41 #include "AliTPCmapper.h"
  42 #include "AliTPCParam.h"
  43 #include "AliTPCCalibRaw.h"
  44 #include "TGraphErrors.h"
  45 #include "AliLog.h"
  46 #include "AliTPCcalibDButil.h"
  47 #include "AliTPCPreprocessorOnline.h"
  48 #include "AliTPCCalibVdrift.h"
  49
  50 ClassImp(AliTPCcalibDButil)
  51 AliTPCcalibDButil::AliTPCcalibDButil() :
  52   TObject(),
  53   fCalibDB(0),
  54   fPadNoise(0x0),
  55   fPedestals(0x0),
  56   fPulserTmean(0x0),
  57   fPulserTrms(0x0),
  58   fPulserQmean(0x0),
  59   fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
  60   fCETmean(0x0),
  61   fCETrms(0x0),
  62   fCEQmean(0x0),
  63   fALTROMasked(0x0),
  64   fCalibRaw(0x0),
  65   fRefPadNoise(0x0),
  66   fRefPedestals(0x0),
  67   fRefPulserTmean(0x0),
  68   fRefPulserTrms(0x0),
  69   fRefPulserQmean(0x0),
  70   fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
  71   fRefCETmean(0x0),
  72   fRefCETrms(0x0),
  73   fRefCEQmean(0x0),
  74   fRefALTROMasked(0x0),
  75   fRefCalibRaw(0x0),
  76   fGoofieArray(0x0),
  77   fMapper(new AliTPCmapper(0x0)),
  78   fNpulserOutliers(-1),
  79   fIrocTimeOffset(0),
  80   fCETmaxLimitAbs(1.5),
  81   fPulTmaxLimitAbs(1.5),
  82   fPulQmaxLimitAbs(5),
  83   fPulQminLimit(11),
  84   fRuns(0),                         // run list with OCDB info
  85   fRunsStart(0),                    // start time for given run
  86   fRunsStop(0)                     // stop time for given run
  87 {
  88   //
  89   // Default ctor
  90   //
  91 }
  92 //_____________________________________________________________________________________
  93 AliTPCcalibDButil::~AliTPCcalibDButil()
  94 {
  95   //
  96   // dtor
  97   //
  98   delete fPulserOutlier;
  99   delete fRefPulserOutlier;
 100   delete fMapper;
 101   if (fRefPadNoise) delete fRefPadNoise;
 102   if (fRefPedestals) delete fRefPedestals;
 103   if (fRefPulserTmean) delete fRefPulserTmean;
 104   if (fRefPulserTrms) delete fRefPulserTrms;
 105   if (fRefPulserQmean) delete fRefPulserQmean;
 106   if (fRefCETmean) delete fRefCETmean;
 107   if (fRefCETrms) delete fRefCETrms;
 108   if (fRefCEQmean) delete fRefCEQmean;
 109   if (fRefALTROMasked) delete fRefALTROMasked;
 110   if (fRefCalibRaw) delete fRefCalibRaw;
 111
 112 }
 113 //_____________________________________________________________________________________
 114 void AliTPCcalibDButil::UpdateFromCalibDB()
 115 {
 116   //
 117   // Update pointers from calibDB
 118   //
 119   if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
 120   fPadNoise=fCalibDB->GetPadNoise();
 121   fPedestals=fCalibDB->GetPedestals();
 122   fPulserTmean=fCalibDB->GetPulserTmean();
 123   fPulserTrms=fCalibDB->GetPulserTrms();
 124   fPulserQmean=fCalibDB->GetPulserQmean();
 125   fCETmean=fCalibDB->GetCETmean();
 126   fCETrms=fCalibDB->GetCETrms();
 127   fCEQmean=fCalibDB->GetCEQmean();
 128   fALTROMasked=fCalibDB->GetALTROMasked();
 129   fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
 130   fCalibRaw=fCalibDB->GetCalibRaw();
 131   UpdatePulserOutlierMap();
 132 }
 133 //_____________________________________________________________________________________
 134 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
 135                                       Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
 136 {
 137   //
 138   // Process the CE data for this run
 139   // the return TVectorD arrays contian the results of the fit
 140   // noutliersCE contains the number of pads marked as outliers,
 141   //   not including masked and edge pads
 142   //
 143
 144   //retrieve CE and ALTRO data
 145   if (!fCETmean){
 146     TString fitString(fitFormula);
 147     fitString.ReplaceAll("++","#");
 148     Int_t ndim=fitString.CountChar('#')+2;
 149     fitResultsA.ResizeTo(ndim);
 150     fitResultsC.ResizeTo(ndim);
 151     fitResultsA.Zero();
 152     fitResultsC.Zero();
 153     noutliersCE=-1;
 154     return;
 155   }
 156   noutliersCE=0;
 157   //create outlier map
 158   AliTPCCalPad *out=0;
 159   if (outCE) out=outCE;
 160   else out=new AliTPCCalPad("outCE","outCE");
 161   AliTPCCalROC *rocMasked=0x0;
 162   //loop over all channels
 163   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 164     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 165     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
 166     AliTPCCalROC *rocOut=out->GetCalROC(iroc);
 167     if (!rocData) {
 168       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
 169       rocOut->Add(1.);
 170       continue;
 171     }
 172     //add time offset to IROCs
 173     if (iroc<AliTPCROC::Instance()->GetNInnerSector())
 174       rocData->Add(fIrocTimeOffset);
 175     //select outliers
 176     UInt_t nrows=rocData->GetNrows();
 177     for (UInt_t irow=0;irow<nrows;++irow){
 178       UInt_t npads=rocData->GetNPads(irow);
 179       for (UInt_t ipad=0;ipad<npads;++ipad){
 180         rocOut->SetValue(irow,ipad,0);
 181         //exclude masked pads
 182         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
 183           rocOut->SetValue(irow,ipad,1);
 184           continue;
 185         }
 186         //exclude first two rows in IROC and last two rows in OROC
 187         if (iroc<36){
 188           if (irow<2) rocOut->SetValue(irow,ipad,1);
 189         } else {
 190           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
 191         }
 192         //exclude edge pads
 193         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
 194         Float_t valTmean=rocData->GetValue(irow,ipad);
 195         //exclude values that are exactly 0
 196         if (valTmean==0) {
 197           rocOut->SetValue(irow,ipad,1);
 198           ++noutliersCE;
 199         }
 200         // exclude channels with too large variations
 201         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
 202           rocOut->SetValue(irow,ipad,1);
 203           ++noutliersCE;
 204         }
 205       }
 206     }
 207   }
 208   //perform fit
 209   TMatrixD dummy;
 210   Float_t chi2Af,chi2Cf;
 211   fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
 212   chi2A=chi2Af;
 213   chi2C=chi2Cf;
 214   if (!outCE) delete out;
 215 }
 216 //_____________________________________________________________________________________
 217 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
 218                      TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
 219                      Float_t &driftTimeA, Float_t &driftTimeC )
 220 {
 221   //
 222   // Calculate statistical information from the CE graphs for drift time and charge
 223   //
 224
 225   //reset arrays
 226   vecTEntries.ResizeTo(72);
 227   vecTMean.ResizeTo(72);
 228   vecTRMS.ResizeTo(72);
 229   vecTMedian.ResizeTo(72);
 230   vecQEntries.ResizeTo(72);
 231   vecQMean.ResizeTo(72);
 232   vecQRMS.ResizeTo(72);
 233   vecQMedian.ResizeTo(72);
 234   vecTEntries.Zero();
 235   vecTMean.Zero();
 236   vecTRMS.Zero();
 237   vecTMedian.Zero();
 238   vecQEntries.Zero();
 239   vecQMean.Zero();
 240   vecQRMS.Zero();
 241   vecQMedian.Zero();
 242   driftTimeA=0;
 243   driftTimeC=0;
 244   TObjArray *arrT=fCalibDB->GetCErocTtime();
 245   TObjArray *arrQ=fCalibDB->GetCErocQtime();
 246   if (arrT){
 247     for (Int_t isec=0;isec<74;++isec){
 248       TGraph *gr=(TGraph*)arrT->At(isec);
 249       if (!gr) continue;
 250       TVectorD values;
 251       Int_t npoints = gr->GetN();
 252       values.ResizeTo(npoints);
 253       Int_t nused =0;
 254       //skip first points, theres always some problems with finding the CE position
 255       for (Int_t ipoint=4; ipoint<npoints; ipoint++){
 256         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
 257           values[nused]=gr->GetY()[ipoint];
 258           nused++;
 259         }
 260       }
 261       //
 262       if (isec<72) vecTEntries[isec]= nused;
 263       if (nused>1){
 264         if (isec<72){
 265           vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 266           vecTMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 267           vecTRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 268         } else if (isec==72){
 269           driftTimeA=TMath::Median(nused,values.GetMatrixArray());
 270         } else if (isec==73){
 271           driftTimeC=TMath::Median(nused,values.GetMatrixArray());
 272         }
 273       }
 274     }
 275   }
 276   if (arrQ){
 277     for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
 278       TGraph *gr=(TGraph*)arrQ->At(isec);
 279       if (!gr) continue;
 280       TVectorD values;
 281       Int_t npoints = gr->GetN();
 282       values.ResizeTo(npoints);
 283       Int_t nused =0;
 284       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 285         if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
 286           values[nused]=gr->GetY()[ipoint];
 287           nused++;
 288         }
 289       }
 290       //
 291       vecQEntries[isec]= nused;
 292       if (nused>1){
 293         vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 294         vecQMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 295         vecQRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 296       }
 297     }
 298   }
 299 }
 300
 301 //_____________________________________________________________________________________
 302 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
 303                       TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
 304                       Int_t &nonMaskedZero)
 305 {
 306   //
 307   // process noise data
 308   // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
 309   //    OROCs small pads [2] and OROCs large pads [3]
 310   // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
 311   // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
 312   //
 313
 314   //set proper size and reset
 315   const UInt_t infoSize=4;
 316   vNoiseMean.ResizeTo(infoSize);
 317   vNoiseMeanSenRegions.ResizeTo(infoSize);
 318   vNoiseRMS.ResizeTo(infoSize);
 319   vNoiseRMSSenRegions.ResizeTo(infoSize);
 320   vNoiseMean.Zero();
 321   vNoiseMeanSenRegions.Zero();
 322   vNoiseRMS.Zero();
 323   vNoiseRMSSenRegions.Zero();
 324   nonMaskedZero=0;
 325   //counters
 326   TVectorD c(infoSize);
 327   TVectorD cs(infoSize);
 328   //tpc parameters
 329   AliTPCParam par;
 330   par.Update();
 331   //retrieve noise and ALTRO data
 332   if (!fPadNoise) return;
 333   AliTPCCalROC *rocMasked=0x0;
 334   //create IROC, OROC1, OROC2 and sensitive region masks
 335   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 336     AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
 337     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 338     UInt_t nrows=noiseROC->GetNrows();
 339     for (UInt_t irow=0;irow<nrows;++irow){
 340       UInt_t npads=noiseROC->GetNPads(irow);
 341       for (UInt_t ipad=0;ipad<npads;++ipad){
 342         //don't use masked channels;
 343         if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
 344         Float_t noiseVal=noiseROC->GetValue(irow,ipad);
 345         //check if noise==0
 346         if (noiseVal==0) {
 347           ++nonMaskedZero;
 348           continue;
 349         }
 350         //check for nan
 351         if ( !(noiseVal<10000000) ){
 352           printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
 353           continue;
 354         }
 355         Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
 356         Int_t masksen=1; // sensitive pards are not masked (0)
 357         if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
 358         if (isec<AliTPCROC::Instance()->GetNInnerSector()){
 359           //IROCs
 360           if (irow>19&&irow<46){
 361             if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
 362           }
 363           Int_t type=1;
 364           vNoiseMean[type]+=noiseVal;
 365           vNoiseRMS[type]+=noiseVal*noiseVal;
 366           ++c[type];
 367           if (!masksen){
 368             vNoiseMeanSenRegions[type]+=noiseVal;
 369             vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 370             ++cs[type];
 371           }
 372         } else {
 373           //OROCs
 374           //define sensive regions
 375           if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
 376           if ( irow>75 ){
 377             Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
 378             if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
 379           }
 380           if ((Int_t)irow<par.GetNRowUp1()){
 381             //OROC1
 382             Int_t type=2;
 383             vNoiseMean[type]+=noiseVal;
 384             vNoiseRMS[type]+=noiseVal*noiseVal;
 385             ++c[type];
 386             if (!masksen){
 387               vNoiseMeanSenRegions[type]+=noiseVal;
 388               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 389               ++cs[type];
 390             }
 391           }else{
 392             //OROC2
 393             Int_t type=3;
 394             vNoiseMean[type]+=noiseVal;
 395             vNoiseRMS[type]+=noiseVal*noiseVal;
 396             ++c[type];
 397             if (!masksen){
 398               vNoiseMeanSenRegions[type]+=noiseVal;
 399               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 400               ++cs[type];
 401             }
 402           }
 403         }
 404         //whole tpc
 405         Int_t type=0;
 406         vNoiseMean[type]+=noiseVal;
 407         vNoiseRMS[type]+=noiseVal*noiseVal;
 408         ++c[type];
 409         if (!masksen){
 410           vNoiseMeanSenRegions[type]+=noiseVal;
 411           vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 412           ++cs[type];
 413         }
 414       }//end loop pads
 415     }//end loop rows
 416   }//end loop sectors (rocs)
 417
 418   //calculate mean and RMS
 419   const Double_t verySmall=0.0000000001;
 420   for (UInt_t i=0;i<infoSize;++i){
 421     Double_t mean=0;
 422     Double_t rms=0;
 423     Double_t meanSen=0;
 424     Double_t rmsSen=0;
 425
 426     if (c[i]>verySmall){
 427 //       printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
 428       mean=vNoiseMean[i]/c[i];
 429       rms=vNoiseRMS[i];
 430       rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
 431     }
 432     vNoiseMean[i]=mean;
 433     vNoiseRMS[i]=rms;
 434
 435     if (cs[i]>verySmall){
 436       meanSen=vNoiseMeanSenRegions[i]/cs[i];
 437       rmsSen=vNoiseRMSSenRegions[i];
 438       rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
 439     }
 440     vNoiseMeanSenRegions[i]=meanSen;
 441     vNoiseRMSSenRegions[i]=rmsSen;
 442   }
 443 }
 444
 445 //_____________________________________________________________________________________
 446 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
 447 {
 448   //
 449   // Process the Pulser information
 450   // vMeanTime:     pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
 451   //
 452
 453   const UInt_t infoSize=4;
 454   //reset counters to error number
 455   vMeanTime.ResizeTo(infoSize);
 456   vMeanTime.Zero();
 457   //counter
 458   TVectorD c(infoSize);
 459   //retrieve pulser and ALTRO data
 460   if (!fPulserTmean) return;
 461   //
 462   //get Outliers
 463   AliTPCCalROC *rocOut=0x0;
 464   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 465     AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
 466     if (!tmeanROC) continue;
 467     rocOut=fPulserOutlier->GetCalROC(isec);
 468     UInt_t nchannels=tmeanROC->GetNchannels();
 469     for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
 470       if (rocOut && rocOut->GetValue(ichannel)) continue;
 471       Float_t val=tmeanROC->GetValue(ichannel);
 472       Int_t type=isec/18;
 473       vMeanTime[type]+=val;
 474       ++c[type];
 475     }
 476   }
 477   //calculate mean
 478   for (UInt_t itype=0; itype<infoSize; ++itype){
 479     if (c[itype]>0) vMeanTime[itype]/=c[itype];
 480     else vMeanTime[itype]=0;
 481   }
 482 }
 483 //_____________________________________________________________________________________
 484 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
 485 {
 486   //
 487   // Get Values from ALTRO configuration data
 488   //
 489   nMasked=-1;
 490   if (!fALTROMasked) return;
 491   nMasked=0;
 492   for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
 493     AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
 494     for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
 495       if (rocMasked->GetValue(ichannel)) ++nMasked;
 496     }
 497   }
 498 }
 499 //_____________________________________________________________________________________
 500 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
 501 {
 502   //
 503   // Proces Goofie values, return statistical information of the currently set goofieArray
 504   // The meaning of the entries are given below
 505   /*
 506   1       TPC_ANODE_I_A00_STAT
 507   2       TPC_DVM_CO2
 508   3       TPC_DVM_DriftVelocity
 509   4       TPC_DVM_FCageHV
 510   5       TPC_DVM_GainFar
 511   6       TPC_DVM_GainNear
 512   7       TPC_DVM_N2
 513   8       TPC_DVM_NumberOfSparks
 514   9       TPC_DVM_PeakAreaFar
 515   10      TPC_DVM_PeakAreaNear
 516   11      TPC_DVM_PeakPosFar
 517   12      TPC_DVM_PeakPosNear
 518   13      TPC_DVM_PickupHV
 519   14      TPC_DVM_Pressure
 520   15      TPC_DVM_T1_Over_P
 521   16      TPC_DVM_T2_Over_P
 522   17      TPC_DVM_T_Over_P
 523   18      TPC_DVM_TemperatureS1
 524    */
 525   if (!fGoofieArray){
 526     Int_t nsensors=19;
 527     vecEntries.ResizeTo(nsensors);
 528     vecMedian.ResizeTo(nsensors);
 529     vecMean.ResizeTo(nsensors);
 530     vecRMS.ResizeTo(nsensors);
 531     vecEntries.Zero();
 532     vecMedian.Zero();
 533     vecMean.Zero();
 534     vecRMS.Zero();
 535     return;
 536   }
 537   Double_t kEpsilon=0.0000000001;
 538   Double_t kBig=100000000000.;
 539   Int_t nsensors = fGoofieArray->NumSensors();
 540   vecEntries.ResizeTo(nsensors);
 541   vecMedian.ResizeTo(nsensors);
 542   vecMean.ResizeTo(nsensors);
 543   vecRMS.ResizeTo(nsensors);
 544   TVectorF values;
 545   for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
 546     AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
 547     if (gsensor &&  gsensor->GetGraph()){
 548       Int_t npoints = gsensor->GetGraph()->GetN();
 549       // filter zeroes
 550       values.ResizeTo(npoints);
 551       Int_t nused =0;
 552       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 553         if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
 554             TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
 555               values[nused]=gsensor->GetGraph()->GetY()[ipoint];
 556               nused++;
 557             }
 558       }
 559       //
 560       vecEntries[isensor]= nused;
 561       if (nused>1){
 562         vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
 563         vecMean[isensor]   = TMath::Mean(nused,values.GetMatrixArray());
 564         vecRMS[isensor]    = TMath::RMS(nused,values.GetMatrixArray());
 565       }
 566     }
 567   }
 568 }
 569 //_____________________________________________________________________________________
 570 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
 571 {
 572   //
 573   // check the variations of the pedestal data to the reference pedestal data
 574   // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
 575   //
 576   const Int_t npar=4;
 577   TVectorF vThres(npar); //thresholds
 578   Int_t nActive=0;       //number of active channels
 579
 580   //reset and set thresholds
 581   pedestalDeviations.ResizeTo(npar);
 582   for (Int_t i=0;i<npar;++i){
 583     pedestalDeviations.GetMatrixArray()[i]=0;
 584     vThres.GetMatrixArray()[i]=(i+1)*.5;
 585   }
 586   //check all needed data is available
 587   if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
 588   //loop over all channels
 589   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 590     AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
 591     AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
 592     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 593     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 594     UInt_t nrows=mROC->GetNrows();
 595     for (UInt_t irow=0;irow<nrows;++irow){
 596       UInt_t npads=mROC->GetNPads(irow);
 597       for (UInt_t ipad=0;ipad<npads;++ipad){
 598         //don't use masked channels;
 599         if (mROC   ->GetValue(irow,ipad)) continue;
 600         if (mRefROC->GetValue(irow,ipad)) continue;
 601         Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
 602         for (Int_t i=0;i<npar;++i){
 603           if (deviation>vThres[i])
 604             ++pedestalDeviations.GetMatrixArray()[i];
 605         }
 606         ++nActive;
 607       }//end ipad
 608     }//ind irow
 609   }//end isec
 610   if (nActive>0){
 611     for (Int_t i=0;i<npar;++i){
 612       pedestalDeviations.GetMatrixArray()[i]/=nActive;
 613     }
 614   }
 615 }
 616 //_____________________________________________________________________________________
 617 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
 618 {
 619   //
 620   // check the variations of the noise data to the reference noise data
 621   // thresholds are 5, 10, 15 and 20 percent respectively.
 622   //
 623   const Int_t npar=4;
 624   TVectorF vThres(npar); //thresholds
 625   Int_t nActive=0;       //number of active channels
 626
 627   //reset and set thresholds
 628   noiseDeviations.ResizeTo(npar);
 629   for (Int_t i=0;i<npar;++i){
 630     noiseDeviations.GetMatrixArray()[i]=0;
 631     vThres.GetMatrixArray()[i]=(i+1)*.05;
 632   }
 633   //check all needed data is available
 634   if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
 635   //loop over all channels
 636   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 637     AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
 638     AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
 639     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 640     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 641     UInt_t nrows=mROC->GetNrows();
 642     for (UInt_t irow=0;irow<nrows;++irow){
 643       UInt_t npads=mROC->GetNPads(irow);
 644       for (UInt_t ipad=0;ipad<npads;++ipad){
 645         //don't use masked channels;
 646         if (mROC   ->GetValue(irow,ipad)) continue;
 647         if (mRefROC->GetValue(irow,ipad)) continue;
 648         Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
 649         for (Int_t i=0;i<npar;++i){
 650           if (deviation>vThres[i])
 651             ++noiseDeviations.GetMatrixArray()[i];
 652         }
 653         ++nActive;
 654       }//end ipad
 655     }//ind irow
 656   }//end isec
 657   if (nActive>0){
 658     for (Int_t i=0;i<npar;++i){
 659       noiseDeviations.GetMatrixArray()[i]/=nActive;
 660     }
 661   }
 662 }
 663 //_____________________________________________________________________________________
 664 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
 665                                                 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
 666 {
 667   //
 668   // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
 669   // thresholds are .5, 1, 5 and 10 percent respectively.
 670   //
 671   //
 672   const Int_t npar=4;
 673   TVectorF vThres(npar); //thresholds
 674   Int_t nActive=0;       //number of active channels
 675
 676   //reset and set thresholds
 677   pulserQdeviations.ResizeTo(npar);
 678   for (Int_t i=0;i<npar;++i){
 679     pulserQdeviations.GetMatrixArray()[i]=0;
 680   }
 681   npadsOutOneTB=0;
 682   npadsOffAdd=0;
 683   varQMean=0;
 684   vThres.GetMatrixArray()[0]=.005;
 685   vThres.GetMatrixArray()[1]=.01;
 686   vThres.GetMatrixArray()[2]=.05;
 687   vThres.GetMatrixArray()[3]=.1;
 688   //check all needed data is available
 689   if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
 690   //
 691   UpdateRefPulserOutlierMap();
 692   //loop over all channels
 693   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 694     AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
 695     AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
 696     AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
 697 //     AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
 698     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 699     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 700     AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
 701     Float_t pt_mean=ptROC->GetMean(oROC);
 702     UInt_t nrows=mROC->GetNrows();
 703     for (UInt_t irow=0;irow<nrows;++irow){
 704       UInt_t npads=mROC->GetNPads(irow);
 705       for (UInt_t ipad=0;ipad<npads;++ipad){
 706         //don't use masked channels;
 707         if (mROC   ->GetValue(irow,ipad)) continue;
 708         if (mRefROC->GetValue(irow,ipad)) continue;
 709         //don't user edge pads
 710         if (ipad==0||ipad==npads-1) continue;
 711         //data
 712         Float_t pq=pqROC->GetValue(irow,ipad);
 713         Float_t pqRef=pqRefROC->GetValue(irow,ipad);
 714         Float_t pt=ptROC->GetValue(irow,ipad);
 715 //         Float_t ptRef=ptRefROC->GetValue(irow,ipad);
 716         //comparisons q
 717         Float_t deviation=TMath::Abs(pq/pqRef-1);
 718         for (Int_t i=0;i<npar;++i){
 719           if (deviation>vThres[i])
 720             ++pulserQdeviations.GetMatrixArray()[i];
 721         }
 722         if (pqRef>11&&pq<11) ++npadsOffAdd;
 723         varQMean+=pq-pqRef;
 724         //comparisons t
 725         if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
 726         ++nActive;
 727       }//end ipad
 728     }//ind irow
 729   }//end isec
 730   if (nActive>0){
 731     for (Int_t i=0;i<npar;++i){
 732       pulserQdeviations.GetMatrixArray()[i]/=nActive;
 733       varQMean/=nActive;
 734     }
 735   }
 736 }
 737 //_____________________________________________________________________________________
 738 void AliTPCcalibDButil::UpdatePulserOutlierMap()
 739 {
 740   //
 741   //
 742   //
 743   PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
 744 }
 745 //_____________________________________________________________________________________
 746 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
 747 {
 748   //
 749   //
 750   //
 751   PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
 752 }
 753 //_____________________________________________________________________________________
 754 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
 755 {
 756   //
 757   // Create a map that contains outliers from the Pulser calibration data.
 758   // The outliers include masked channels, edge pads and pads with
 759   //   too large timing and charge variations.
 760   // fNpulserOutliers is the number of outliers in the Pulser calibration data.
 761   //   those do not contain masked and edge pads
 762   //
 763   if (!pulT||!pulQ) {
 764     //reset map
 765     pulOut->Multiply(0.);
 766     fNpulserOutliers=-1;
 767     return;
 768   }
 769   AliTPCCalROC *rocMasked=0x0;
 770   fNpulserOutliers=0;
 771
 772   //Create Outlier Map
 773   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 774     AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
 775     AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
 776     AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
 777     if (!tmeanROC||!qmeanROC) {
 778       //reset outliers in this ROC
 779       outROC->Multiply(0.);
 780       continue;
 781     }
 782     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 783 //     Double_t dummy=0;
 784 //     Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
 785 //     Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
 786     UInt_t nrows=tmeanROC->GetNrows();
 787     for (UInt_t irow=0;irow<nrows;++irow){
 788       UInt_t npads=tmeanROC->GetNPads(irow);
 789       for (UInt_t ipad=0;ipad<npads;++ipad){
 790         Int_t outlier=0,masked=0;
 791         Float_t q=qmeanROC->GetValue(irow,ipad);
 792         Float_t t=tmeanROC->GetValue(irow,ipad);
 793         //masked channels are outliers
 794         if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
 795         //edge pads are outliers
 796         if (ipad==0||ipad==npads-1) masked=1;
 797         //channels with too large charge or timing deviation from the meadian are outliers
 798 //         if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
 799         if (q<fPulQminLimit && !masked) outlier=1;
 800         //check for nan
 801         if ( !(q<10000000) || !(t<10000000)) outlier=1;
 802         outROC->SetValue(irow,ipad,outlier+masked);
 803         fNpulserOutliers+=outlier;
 804       }
 805     }
 806   }
 807 }
 808 //_____________________________________________________________________________________
 809 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
 810 {
 811   //
 812   // Create pad time0 object from pulser and/or CE data, depending on the selected model
 813   // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
 814   // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
 815   // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
 816   //
 817   // In case model 2 is invoked - gy arival time gradient is also returned
 818   //
 819   gyA=0;
 820   gyC=0;
 821   AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
 822   // decide between different models
 823   if (model==0||model==1){
 824     TVectorD vMean;
 825     if (model==1) ProcessPulser(vMean);
 826     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 827       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 828       if (!rocPulTmean) continue;
 829       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 830       AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
 831       Float_t mean=rocPulTmean->GetMean(rocOut);
 832       //treat case where a whole partition is masked
 833       if (mean==0) mean=rocPulTmean->GetMean();
 834       if (model==1) {
 835         Int_t type=isec/18;
 836         mean=vMean[type];
 837       }
 838       UInt_t nrows=rocTime0->GetNrows();
 839       for (UInt_t irow=0;irow<nrows;++irow){
 840         UInt_t npads=rocTime0->GetNPads(irow);
 841         for (UInt_t ipad=0;ipad<npads;++ipad){
 842           Float_t time=rocPulTmean->GetValue(irow,ipad);
 843           //in case of an outlier pad use the mean of the altro values.
 844           //This should be the most precise guess in that case.
 845           if (rocOut->GetValue(irow,ipad)) {
 846             time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
 847             if (time==0) time=mean;
 848           }
 849           Float_t val=time-mean;
 850           rocTime0->SetValue(irow,ipad,val);
 851         }
 852       }
 853     }
 854   } else if (model==2){
 855     Double_t pgya,pgyc,pchi2a,pchi2c;
 856     AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
 857     fCETmean->Add(padPulser,-1.);
 858     TVectorD vA,vC;
 859     AliTPCCalPad outCE("outCE","outCE");
 860     Int_t nOut;
 861     ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
 862     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
 863 //     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
 864     if (!padFit) { delete padPulser; return 0;}
 865     gyA=vA[2];
 866     gyC=vC[2];
 867     fCETmean->Add(padPulser,1.);
 868     padTime0->Add(fCETmean);
 869     padTime0->Add(padFit,-1);
 870     delete padPulser;
 871     TVectorD vFitROC;
 872     TMatrixD mFitROC;
 873     Float_t chi2;
 874     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 875       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 876       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 877       AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
 878       AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
 879       rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
 880       AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
 881       Float_t mean=rocPulTmean->GetMean(rocOutPul);
 882       if (mean==0) mean=rocPulTmean->GetMean();
 883       UInt_t nrows=rocTime0->GetNrows();
 884       for (UInt_t irow=0;irow<nrows;++irow){
 885         UInt_t npads=rocTime0->GetNPads(irow);
 886         for (UInt_t ipad=0;ipad<npads;++ipad){
 887           Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
 888           if (rocOutCE->GetValue(irow,ipad)){
 889             Float_t valOut=rocCEfit->GetValue(irow,ipad);
 890             if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
 891             rocTime0->SetValue(irow,ipad,valOut);
 892           }
 893         }
 894       }
 895       delete rocCEfit;
 896     }
 897     delete padFit;
 898   }
 899   Double_t median = padTime0->GetMedian();
 900   padTime0->Add(-median);  // normalize to median
 901   return padTime0;
 902 }
 903 //_____________________________________________________________________________________
 904 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
 905 {
 906   if (roc==0) return 0.;
 907   const Int_t sector=roc->GetSector();
 908   AliTPCROC *tpcRoc=AliTPCROC::Instance();
 909   const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
 910   Float_t mean=0;
 911   Int_t   n=0;
 912
 913   //loop over a small range around the requested pad (+-10 rows/pads)
 914   for (Int_t irow=row-10;irow<row+10;++irow){
 915     if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
 916     for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
 917       if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
 918       const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
 919       if (altroRoc!=altroCurr) continue;
 920       if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
 921       Float_t val=roc->GetValue(irow,ipad);
 922       mean+=val;
 923       ++n;
 924     }
 925   }
 926   if (n>0) mean/=n;
 927   return mean;
 928 }
 929 //_____________________________________________________________________________________
 930 void AliTPCcalibDButil::SetRefFile(const char* filename)
 931 {
 932   //
 933   // load cal pad objects form the reference file
 934   //
 935   TDirectory *currDir=gDirectory;
 936   TFile f(filename);
 937   fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
 938   fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
 939   //pulser data
 940   fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
 941   fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
 942   fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
 943   //CE data
 944   fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
 945   fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
 946   fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
 947   //Altro data
 948 //   fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
 949 //   fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
 950 //   fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
 951 //   fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
 952   fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
 953   f.Close();
 954   currDir->cd();
 955 }
 956
 957
 958
 959
 960 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin,  Float_t cutTrmsMax, Float_t cutMaxDistT){
 961   //
 962   // Author:  marian.ivanov@cern.ch
 963   //
 964   // Create outlier map for CE study
 965   // Parameters:
 966   //  Return value - outlyer map
 967   //  noutlyersCE  - number of outlyers
 968   //  minSignal    - minimal total Q signal
 969   //  cutRMSMin    - minimal width of the signal in respect to the median
 970   //  cutRMSMax    - maximal width of the signal in respect to the median
 971   //  cutMaxDistT  - maximal deviation from time median per chamber
 972   //
 973   // Outlyers criteria:
 974   // 0. Exclude masked pads
 975   // 1. Exclude first two rows in IROC and last two rows in OROC
 976   // 2. Exclude edge pads
 977   // 3. Exclude channels with too large variations
 978   // 4. Exclude pads with too small signal
 979   // 5. Exclude signal with outlyers RMS
 980   // 6. Exclude channels to far from the chamber median
 981   noutliersCE=0;
 982   //create outlier map
 983   AliTPCCalPad *out=ceOut;
 984   if (!out)     out= new AliTPCCalPad("outCE","outCE");
 985   AliTPCCalROC *rocMasked=0x0;
 986   if (!fCETmean) return 0;
 987   if (!fCETrms) return 0;
 988   if (!fCEQmean) return 0;
 989   //
 990   //loop over all channels
 991   //
 992   Double_t rmsMedian         = fCETrms->GetMedian();
 993   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 994     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 995     if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
 996     AliTPCCalROC *rocOut       = out->GetCalROC(iroc);
 997     AliTPCCalROC *rocCEQ       = fCEQmean->GetCalROC(iroc);
 998     AliTPCCalROC *rocCETrms    = fCETrms->GetCalROC(iroc);
 999     Double_t trocMedian        = rocData->GetMedian();
1000     //
1001     if (!rocData) {
1002       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1003       rocOut->Add(1.);
1004       continue;
1005     }
1006     //
1007     //select outliers
1008     UInt_t nrows=rocData->GetNrows();
1009     for (UInt_t irow=0;irow<nrows;++irow){
1010       UInt_t npads=rocData->GetNPads(irow);
1011       for (UInt_t ipad=0;ipad<npads;++ipad){
1012         rocOut->SetValue(irow,ipad,0);
1013         Float_t valTmean=rocData->GetValue(irow,ipad);
1014         Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1015         Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1016         //0. exclude masked pads
1017         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1018           rocOut->SetValue(irow,ipad,1);
1019           continue;
1020         }
1021         //1. exclude first two rows in IROC and last two rows in OROC
1022         if (iroc<36){
1023           if (irow<2) rocOut->SetValue(irow,ipad,1);
1024         } else {
1025           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1026         }
1027         //2. exclude edge pads
1028         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1029         //exclude values that are exactly 0
1030         if (valTmean==0) {
1031           rocOut->SetValue(irow,ipad,1);
1032           ++noutliersCE;
1033         }
1034         //3.  exclude channels with too large variations
1035         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1036           rocOut->SetValue(irow,ipad,1);
1037           ++noutliersCE;
1038         }
1039         //
1040         //4.  exclude channels with too small signal
1041         if (valQmean<minSignal) {
1042           rocOut->SetValue(irow,ipad,1);
1043           ++noutliersCE;
1044         }
1045         //
1046         //5. exclude channels with too small rms
1047         if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1048           rocOut->SetValue(irow,ipad,1);
1049           ++noutliersCE;
1050         }
1051         //
1052         //6. exclude channels to far from the chamber median
1053         if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1054           rocOut->SetValue(irow,ipad,1);
1055           ++noutliersCE;
1056         }
1057       }
1058     }
1059   }
1060   //
1061   return out;
1062 }
1063
1064
1065 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1066   //
1067   // Author: marian.ivanov@cern.ch
1068   //
1069   // Create outlier map for Pulser
1070   // Parameters:
1071   //  Return value     - outlyer map
1072   //  noutlyersPulser  - number of outlyers
1073   //  cutTime          - absolute cut - distance to the median of chamber
1074   //  cutnRMSQ         - nsigma cut from median  q distribution per chamber
1075   //  cutnRMSrms       - nsigma cut from median  rms distribution
1076   // Outlyers criteria:
1077   // 0. Exclude masked pads
1078   // 1. Exclude time outlyers (default 3 time bins)
1079   // 2. Exclude q outlyers    (default 5 sigma)
1080   // 3. Exclude rms outlyers  (default 5 sigma)
1081   noutliersPulser=0;
1082   AliTPCCalPad *out=pulserOut;
1083   if (!out)     out= new AliTPCCalPad("outPulser","outPulser");
1084   AliTPCCalROC *rocMasked=0x0;
1085   if (!fPulserTmean) return 0;
1086   if (!fPulserTrms) return 0;
1087   if (!fPulserQmean) return 0;
1088   //
1089   //loop over all channels
1090   //
1091   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1092     if (fALTROMasked)   rocMasked= fALTROMasked->GetCalROC(iroc);
1093     AliTPCCalROC *rocData       = fPulserTmean->GetCalROC(iroc);
1094     AliTPCCalROC *rocOut        = out->GetCalROC(iroc);
1095     AliTPCCalROC *rocPulserQ    = fPulserQmean->GetCalROC(iroc);
1096     AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1097     //
1098     Double_t rocMedianT         = rocData->GetMedian();
1099     Double_t rocMedianQ         = rocPulserQ->GetMedian();
1100     Double_t rocRMSQ            = rocPulserQ->GetRMS();
1101     Double_t rocMedianTrms      = rocPulserTrms->GetMedian();
1102     Double_t rocRMSTrms         = rocPulserTrms->GetRMS();
1103     for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1104       rocOut->SetValue(ichannel,0);
1105       Float_t valTmean=rocData->GetValue(ichannel);
1106       Float_t valQmean=rocPulserQ->GetValue(ichannel);
1107       Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1108       Int_t isOut=0;
1109       if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1110       if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1111       if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1112       rocOut->SetValue(ichannel,isOut);
1113       if (isOut) noutliersPulser++;
1114     }
1115   }
1116   return out;
1117 }
1118
1119
1120 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1121   //
1122   // Author : Marian Ivanov
1123   // Create pad time0 correction map using information from the CE and from pulser
1124   //
1125   //
1126   // Return PadTime0 to be used for time0 relative alignment
1127   // if dump file specified intermediat results are dumped to the fiel and can be visualized
1128   // using $ALICE_ROOT/TPC/script/gui application
1129   //
1130   // fitResultsA - fitParameters A side
1131   // fitResultsC - fitParameters C side
1132   // chi2A       - chi2/ndf for A side (assuming error 1 time bin)
1133   // chi2C       - chi2/ndf for C side (assuming error 1 time bin)
1134   //
1135   //
1136   // Algorithm:
1137   // 1. Find outlier map for CE
1138   // 2. Find outlier map for Pulser
1139   // 3. Replace outlier by median at given sector  (median without outliers)
1140   // 4. Substract from the CE data pulser
1141   // 5. Fit the CE with formula
1142   //    5.1) (IROC-OROC) offset
1143   //    5.2) gx
1144   //    5.3) gy
1145   //    5.4) (lx-xmid)
1146   //    5.5) (IROC-OROC)*(lx-xmid)
1147   //    5.6) (ly/lx)^2
1148   // 6. Substract gy fit dependence from the CE data
1149   // 7. Add pulser back to CE data
1150   // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1151   // 9. return CE data
1152   //
1153   // Time0 <= padCE = padCEin  -padCEfitGy  - if not outlier
1154   // Time0 <= padCE = padFitAll-padCEfitGy  - if outlier
1155
1156   // fit formula
1157   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)";
1158   // output for fit formula
1159   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)";
1160   // gy part of formula
1161   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)";
1162   //
1163   //
1164   if (!fCETmean) return 0;
1165   Double_t pgya,pgyc,pchi2a,pchi2c;
1166   AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1167   AliTPCCalPad * padCEOut     = CreateCEOutlyerMap(nOut);
1168
1169   AliTPCCalPad * padPulser    = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1170   AliTPCCalPad * padCE        = new AliTPCCalPad(*fCETmean);
1171   AliTPCCalPad * padCEIn      = new AliTPCCalPad(*fCETmean);
1172   AliTPCCalPad * padOut       = new AliTPCCalPad("padOut","padOut");
1173   padPulser->SetName("padPulser");
1174   padPulserOut->SetName("padPulserOut");
1175   padCE->SetName("padCE");
1176   padCEIn->SetName("padCEIn");
1177   padCEOut->SetName("padCEOut");
1178   padOut->SetName("padOut");
1179
1180   //
1181   // make combined outlyers map
1182   // and replace outlyers in maps with median for chamber
1183   //
1184   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1185     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1186     AliTPCCalROC * rocPulser    = padPulser->GetCalROC(iroc);
1187     AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1188     AliTPCCalROC * rocCEOut     = padCEOut->GetCalROC(iroc);
1189     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1190     Double_t ceMedian           = rocCE->GetMedian(rocCEOut);
1191     Double_t pulserMedian       = rocPulser->GetMedian(rocCEOut);
1192     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1193       if (rocPulserOut->GetValue(ichannel)>0) {
1194         rocPulser->SetValue(ichannel,pulserMedian);
1195         rocOut->SetValue(ichannel,1);
1196       }
1197       if (rocCEOut->GetValue(ichannel)>0) {
1198         rocCE->SetValue(ichannel,ceMedian);
1199         rocOut->SetValue(ichannel,1);
1200       }
1201     }
1202   }
1203   //
1204   // remove pulser time 0
1205   //
1206   padCE->Add(padPulser,-1);
1207   //
1208   // Make fits
1209   //
1210   TMatrixD dummy;
1211   Float_t chi2Af,chi2Cf;
1212   padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1213   chi2A=chi2Af;
1214   chi2C=chi2Cf;
1215   //
1216   AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1217   padCEFitGY->SetName("padCEFitGy");
1218   //
1219   AliTPCCalPad *padCEFit  =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1220   padCEFit->SetName("padCEFit");
1221   //
1222   AliTPCCalPad* padCEDiff  = new AliTPCCalPad(*padCE);
1223   padCEDiff->SetName("padCEDiff");
1224   padCEDiff->Add(padCEFit,-1.);
1225   //
1226   //
1227   padCE->Add(padCEFitGY,-1.);
1228
1229   padCE->Add(padPulser,1.);
1230   Double_t padmedian = padCE->GetMedian();
1231   padCE->Add(-padmedian);  // normalize to median
1232   //
1233   // Replace outliers by fit value - median of diff per given chamber -GY fit
1234   //
1235   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1236     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1237     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1238     AliTPCCalROC * rocCEFit     = padCEFit->GetCalROC(iroc);
1239     AliTPCCalROC * rocCEFitGY   = padCEFitGY->GetCalROC(iroc);
1240     AliTPCCalROC * rocCEDiff    = padCEDiff->GetCalROC(iroc);
1241     //
1242     Double_t diffMedian         = rocCEDiff->GetMedian(rocOut);
1243     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1244       if (rocOut->GetValue(ichannel)==0) continue;
1245       Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1246       rocCE->SetValue(ichannel,value);
1247     }
1248   }
1249   //
1250   //
1251   if (dumpfile){
1252     //dump to the file - result can be visualized
1253     AliTPCPreprocessorOnline preprocesor;
1254     preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1255     preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1256     preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1257     preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1258     //
1259     preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1260     preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1261     //
1262     preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1263     preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1264     preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1265     preprocesor.DumpToFile(dumpfile);
1266   }
1267   delete padPulser;
1268   delete padPulserOut;
1269   delete padCEIn;
1270   delete padCEOut;
1271   delete padOut;
1272   delete padCEDiff;
1273   delete padCEFitGY;
1274   return padCE;
1275 }
1276
1277
1278
1279
1280
1281 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1282   //
1283   // find the closest point to xref  in x  direction
1284   // return dx and value
1285   Int_t index=0;
1286   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1287   if (index<0) index=0;
1288   if (index>=graph->GetN()-1) index=graph->GetN()-2;
1289   if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1290   dx = xref-graph->GetX()[index];
1291   y  = graph->GetY()[index];
1292   return index;
1293 }
1294
1295
1296 Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1297   //
1298   // Get the correction of the trigger offset
1299   // combining information from the laser track calibration
1300   // and from cosmic calibration
1301   //
1302   // run       - run number
1303   // timeStamp - tim stamp in seconds
1304   // deltaT    - integration period to calculate offset
1305   // deltaTLaser -max validity of laser data
1306   // valType   - 0 - median, 1- mean
1307   //
1308   // Integration vaues are just recomendation - if not possible to get points
1309   // automatically increase the validity by factor 2
1310   // (recursive algorithm until one month of data taking)
1311   //
1312   //
1313   const Float_t kLaserCut=0.0005;
1314   const Int_t   kMaxPeriod=3600*24*30*3; // 3 month max
1315   const Int_t   kMinPoints=20;
1316   //
1317   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1318   if (!array) {
1319     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1320   }
1321   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1322   if (!array) return 0;
1323   //
1324   TGraphErrors *laserA[3]={0,0,0};
1325   TGraphErrors *laserC[3]={0,0,0};
1326   TGraphErrors *cosmicAll=0;
1327   laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1328   laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1329   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1330   //
1331   //
1332   if (!cosmicAll) return 0;
1333   Int_t nmeasC=cosmicAll->GetN();
1334   Float_t *tdelta = new Float_t[nmeasC];
1335   Int_t nused=0;
1336   for (Int_t i=0;i<nmeasC;i++){
1337     if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1338     Float_t ccosmic=cosmicAll->GetY()[i];
1339     Double_t yA=0,yC=0,dA=0,dC=0;
1340     if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1341     if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1342     //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1343     //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1344     //
1345     if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1346     Float_t claser=0;
1347     if (TMath::Abs(yA-yC)<kLaserCut) {
1348       claser=(yA-yC)*0.5;
1349     }else{
1350       if (i%2==0)  claser=yA;
1351       if (i%2==1)  claser=yC;
1352     }
1353     tdelta[nused]=ccosmic-claser;
1354     nused++;
1355   }
1356   if (nused<kMinPoints &&deltaT<kMaxPeriod) return  AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1357   Double_t median = TMath::Median(nused,tdelta);
1358   Double_t mean  = TMath::Mean(nused,tdelta);
1359   delete tdelta;
1360   return (valType==0) ? median:mean;
1361 }
1362
1363 Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1364   //
1365   // Get the correction of the drift velocity
1366   // combining information from the laser track calibration
1367   // and from cosmic calibration
1368   //
1369   // dist      - return value - distance to closest point in graph
1370   // run       - run number
1371   // timeStamp - tim stamp in seconds
1372   // deltaT    - integration period to calculate time0 offset
1373   // deltaTLaser -max validity of laser data
1374   // valType   - 0 - median, 1- mean
1375   //
1376   // Integration vaues are just recomendation - if not possible to get points
1377   // automatically increase the validity by factor 2
1378   // (recursive algorithm until one month of data taking)
1379   //
1380   //
1381   //
1382   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1383   if (!array) {
1384     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1385   }
1386   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1387   if (!array) return 0;
1388   TGraphErrors *cosmicAll=0;
1389   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1390   if (!cosmicAll) return 0;
1391   Double_t grY=0;
1392   AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1393
1394   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1395   Double_t vcosmic=  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1396   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1397   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
1398   return  vcosmic+t0;
1399
1400   /*
1401     Example usage:
1402
1403     Int_t run=89000
1404     TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1405     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1406     laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1407     //
1408     Double_t *yvd= new Double_t[cosmicAll->GetN()];
1409     Double_t *yt0= new Double_t[cosmicAll->GetN()];
1410     for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1411     for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1412
1413     TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1414     TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1415
1416   */
1417
1418 }
1419
1420 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1421   //
1422   // Get the correction of the drift velocity using the laser tracks calbration
1423   //
1424   // run       - run number
1425   // timeStamp - tim stamp in seconds
1426   // deltaT    - integration period to calculate time0 offset
1427   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1428   // Note in case no data form both A and C side - the value from active side used
1429   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1430   TGraphErrors *grlaserA=0;
1431   TGraphErrors *grlaserC=0;
1432   Double_t vlaserA=0, vlaserC=0;
1433   if (!array) return 0;
1434   grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1435   grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1436   Double_t deltaY;
1437   if (grlaserA) {
1438     AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1439     if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
1440     else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1441   }
1442   if (grlaserC) {
1443     AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1444     if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
1445     else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1446   }
1447   if (side==0) return vlaserA;
1448   if (side==1) return vlaserC;
1449   Double_t mdrift=(vlaserA+vlaserC)*0.5;
1450   if (!grlaserA) return vlaserC;
1451   if (!grlaserC) return vlaserA;
1452   return mdrift;
1453 }
1454
1455
1456 Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1457   //
1458   // Get the correction of the drift velocity using the CE laser data
1459   // combining information from the CE,  laser track calibration
1460   // and P/T calibration
1461   //
1462   // run       - run number
1463   // timeStamp - tim stamp in seconds
1464   // deltaT    - integration period to calculate time0 offset
1465   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1466   TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
1467   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
1468   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1469   AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1470   //
1471   //
1472   Double_t corrPTA = 0, corrPTC=0;
1473   Double_t ltime0A = 0, ltime0C=0;
1474   Double_t gry=0;
1475   Double_t corrA=0, corrC=0;
1476   Double_t timeA=0, timeC=0;
1477   TGraph *graphA = (TGraph*)arrT->At(72);
1478   TGraph *graphC = (TGraph*)arrT->At(73);
1479   if (!graphA && !graphC) return 0.;
1480   if (graphA &&graphA->GetN()>0) {
1481     AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1482     timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1483     Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1484     ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1485     if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
1486     corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1487     corrA-=corrPTA;
1488   }
1489   if (graphC&&graphC->GetN()>0){
1490     AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1491     timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1492     Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1493     ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1494     if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
1495     corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1496     corrC-=corrPTC;
1497   }
1498
1499   if (side ==0 ) return corrA;
1500   if (side ==1 ) return corrC;
1501   Double_t corrM= (corrA+corrC)*0.5;
1502   if (!graphA) corrM=corrC;
1503   if (!graphC) corrM=corrA;
1504   return corrM;
1505 }
1506
1507
1508
1509
1510 Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1511   //
1512   // VERY obscure method - we need something in framework
1513   // Find the TPC runs with temperature OCDB entry
1514   // cache the start and end of the run
1515   //
1516   AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1517   if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1518   if (!storage) return 0;
1519   TString path=storage->GetURI();
1520   TString runsT;
1521   {
1522     TString command;
1523     if (path.Contains("local")){  // find the list if local system
1524       path.ReplaceAll("local://","");
1525       path+="TPC/Calib/Temperature";
1526       command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
1527     }
1528     runsT=gSystem->GetFromPipe(command);
1529   }
1530   TObjArray *arr= runsT.Tokenize("r");
1531   if (!arr) return 0;
1532   //
1533   TArrayI indexes(arr->GetEntries());
1534   TArrayI runs(arr->GetEntries());
1535   Int_t naccept=0;
1536   {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1537       Int_t irunN = atoi(arr->At(irun)->GetName());
1538       if (irunN<startRun) continue;
1539       if (irunN>stopRun) continue;
1540       runs[naccept]=irunN;
1541       naccept++;
1542     }}
1543   fRuns.Set(naccept);
1544   fRunsStart.Set(fRuns.fN);
1545   fRunsStop.Set(fRuns.fN);
1546   TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1547   for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
1548
1549   //
1550   AliCDBEntry * entry = 0;
1551   {for (Int_t irun=0;irun<fRuns.fN; irun++){
1552       entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1553       if (!entry) continue;
1554       AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1555       if (!tmpRun) continue;
1556       fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1557       fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1558       //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1559     }}
1560   return fRuns.fN;
1561 }
1562
1563
1564 Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
1565   //
1566   // binary search - find the run for given time stamp
1567   //
1568   Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
1569   Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
1570   Int_t cindex  = -1;
1571   for (Int_t index=index0; index<=index1; index++){
1572     if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
1573     if (debug) {
1574       printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
1575     }
1576   }
1577   if (cindex<0) cindex =(index0+index1)/2;
1578   if (cindex<0) {
1579     return 0;
1580   }
1581   return fRuns[cindex];
1582 }
1583
1584
1585
1586
1587
1588 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
1589   //
1590   // filter outlyer measurement
1591   // Only points around median +- sigmaCut filtered
1592   if (!graph) return  0;
1593   Int_t kMinPoints=2;
1594   Int_t npoints0 = graph->GetN();
1595   Int_t npoints=0;
1596   Float_t  rmsY=0;
1597   Double_t *outx=new Double_t[npoints0];
1598   Double_t *outy=new Double_t[npoints0];
1599   //
1600   //
1601   if (npoints0<kMinPoints) return 0;
1602   for (Int_t iter=0; iter<3; iter++){
1603     npoints=0;
1604     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1605       if (graph->GetY()[ipoint]==0) continue;
1606       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
1607       outx[npoints]  = graph->GetX()[ipoint];
1608       outy[npoints]  = graph->GetY()[ipoint];
1609       npoints++;
1610     }
1611     if (npoints<=1) break;
1612     medianY  =TMath::Median(npoints,outy);
1613     rmsY   =TMath::RMS(npoints,outy);
1614   }
1615   TGraph *graphOut=0;
1616   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1617   return graphOut;
1618 }
1619
1620
1621 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
1622   //
1623   // filter outlyer measurement
1624   // Only points around median +- cut filtered
1625   if (!graph) return  0;
1626   Int_t kMinPoints=2;
1627   Int_t npoints0 = graph->GetN();
1628   Int_t npoints=0;
1629   Float_t  rmsY=0;
1630   Double_t *outx=new Double_t[npoints0];
1631   Double_t *outy=new Double_t[npoints0];
1632   //
1633   //
1634   if (npoints0<kMinPoints) return 0;
1635   for (Int_t iter=0; iter<3; iter++){
1636     npoints=0;
1637     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1638       if (graph->GetY()[ipoint]==0) continue;
1639       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
1640       outx[npoints]  = graph->GetX()[ipoint];
1641       outy[npoints]  = graph->GetY()[ipoint];
1642       npoints++;
1643     }
1644     if (npoints<=1) break;
1645     medianY  =TMath::Median(npoints,outy);
1646     rmsY   =TMath::RMS(npoints,outy);
1647   }
1648   TGraph *graphOut=0;
1649   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1650   return graphOut;
1651 }
1652
1653
1654
1655 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
1656   //
1657   // filter outlyer measurement
1658   // Only points with normalized errors median +- sigmaCut filtered
1659   //
1660   Int_t kMinPoints=10;
1661   Int_t npoints0 = graph->GetN();
1662   Int_t npoints=0;
1663   Float_t  medianErr=0, rmsErr=0;
1664   Double_t *outx=new Double_t[npoints0];
1665   Double_t *outy=new Double_t[npoints0];
1666   Double_t *erry=new Double_t[npoints0];
1667   Double_t *nerry=new Double_t[npoints0];
1668   Double_t *errx=new Double_t[npoints0];
1669   //
1670   //
1671   if (npoints0<kMinPoints) return 0;
1672   for (Int_t iter=0; iter<3; iter++){
1673     npoints=0;
1674     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1675       nerry[npoints]  = graph->GetErrorY(ipoint);
1676       if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
1677       erry[npoints]  = graph->GetErrorY(ipoint);
1678       outx[npoints]  = graph->GetX()[ipoint];
1679       outy[npoints]  = graph->GetY()[ipoint];
1680       errx[npoints]  = graph->GetErrorY(ipoint);
1681       npoints++;
1682     }
1683     if (npoints==0) break;
1684     medianErr=TMath::Median(npoints,erry);
1685     medianY  =TMath::Median(npoints,outy);
1686     rmsErr   =TMath::RMS(npoints,erry);
1687   }
1688   TGraphErrors *graphOut=0;
1689   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
1690   delete []outx;
1691   delete []outy;
1692   delete []errx;
1693   delete []erry;
1694   return graphOut;
1695 }
1696
1697
1698 void AliTPCcalibDButil::Sort(TGraph *graph){
1699   //
1700   // sort array - neccessay for approx
1701   //
1702   Int_t npoints = graph->GetN();
1703   Int_t *indexes=new Int_t[npoints];
1704   Double_t *outx=new Double_t[npoints];
1705   Double_t *outy=new Double_t[npoints];
1706   TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
1707   for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
1708   for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
1709   for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
1710   for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
1711
1712 }
1713 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
1714   //
1715   // smmoth graph - mean on the interval
1716   //
1717   Sort(graph);
1718   Int_t npoints = graph->GetN();
1719   Double_t *outy=new Double_t[npoints];
1720   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1721     Double_t lx=graph->GetX()[ipoint];
1722     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1723     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1724     if (index0<0) index0=0;
1725     if (index1>=npoints-1) index1=npoints-1;
1726     if ((index1-index0)>1){
1727       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1728     }else{
1729       outy[ipoint]=graph->GetY()[ipoint];
1730     }
1731   }
1732   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1733     graph->GetY()[ipoint] = outy[ipoint];
1734   }
1735   delete[] outy;
1736 }
1737
1738 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
1739   //
1740   // Use constant interpolation outside of range
1741   //
1742   if (!graph) {
1743     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
1744     return 0;
1745   }
1746   if (graph->GetN()<1){
1747     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
1748     return 0;
1749   }
1750   if (xref<graph->GetX()[0]) return graph->GetY()[0];
1751   if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
1752   return graph->Eval( xref);
1753 }
1754
1755 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy,  Double_t sigmaCut){
1756   //
1757   // Filter DCS sensor information
1758   //   ymin     - minimal value
1759   //   ymax     - max value
1760   //   maxdy    - maximal deirivative
1761   //   sigmaCut - cut on values and derivative in terms of RMS distribution
1762   // Return value - accepted fraction
1763   //
1764   // Algorithm:
1765   //
1766   // 0. Calculate median and rms of values in specified range
1767   // 1. Filter out outliers - median+-sigmaCut*rms
1768   //    values replaced by median
1769   //
1770   AliSplineFit * fit    = sensor->GetFit();
1771   if (!fit) return 0.;
1772   Int_t          nknots = fit->GetKnots();
1773   if (nknots==0) {
1774     delete fit;
1775     sensor->SetFit(0);
1776     return 0;
1777   }
1778   //
1779   Double_t *yin0  = new Double_t[nknots];
1780   Double_t *yin1  = new Double_t[nknots];
1781   Int_t naccept=0;
1782
1783   for (Int_t iknot=0; iknot< nknots; iknot++){
1784     if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
1785       yin0[naccept]  = fit->GetY0()[iknot];
1786       yin1[naccept]  = fit->GetY1()[iknot];
1787       if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
1788       naccept++;
1789     }
1790   }
1791   if (naccept<1) {
1792     delete fit;
1793     sensor->SetFit(0);
1794     return 0.;
1795   }
1796   Double_t medianY0=0, medianY1=0;
1797   Double_t rmsY0   =0, rmsY1=0;
1798   medianY0 = TMath::Median(naccept, yin0);
1799   medianY1 = TMath::Median(naccept, yin1);
1800   rmsY0    = TMath::RMS(naccept, yin0);
1801   rmsY1    = TMath::RMS(naccept, yin1);
1802   naccept=0;
1803   //
1804   // 1. Filter out outliers - median+-sigmaCut*rms
1805   //    values replaced by median
1806   //    if replaced the derivative set to 0
1807   //
1808   for (Int_t iknot=0; iknot< nknots; iknot++){
1809     Bool_t isOK=kTRUE;
1810     if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
1811     if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
1812     if (nknots<2) fit->GetY1()[iknot]=0;
1813     if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
1814     if (!isOK){
1815       fit->GetY0()[iknot]=medianY0;
1816       fit->GetY1()[iknot]=0;
1817     }else{
1818       naccept++;
1819     }
1820   }
1821   delete [] yin0;
1822   delete [] yin1;
1823   return Float_t(naccept)/Float_t(nknots);
1824 }
1825
1826 Float_t  AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
1827   //
1828   // Filter temperature array
1829   // tempArray    - array of temperatures         -
1830   // ymin         - minimal accepted temperature  - default 15
1831   // ymax         - maximal accepted temperature  - default 22
1832   // sigmaCut     - values filtered on interval median+-sigmaCut*rms - defaut 5
1833   // return value - fraction of filtered sensors
1834   const Double_t kMaxDy=0.1;
1835   Int_t nsensors=tempArray->NumSensors();
1836   if (nsensors==0) return 0.;
1837   Int_t naccept=0;
1838   for (Int_t isensor=0; isensor<nsensors; isensor++){
1839     AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
1840     if (!sensor) continue;
1841     //printf("%d\n",isensor);
1842     FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
1843     if (sensor->GetFit()==0){
1844       delete sensor;
1845       tempArray->RemoveSensorNum(isensor);
1846     }else{
1847       naccept++;
1848     }
1849   }
1850   return Float_t(naccept)/Float_t(nsensors);
1851 }
1852
1853
1854 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
1855   //
1856   // Filter CE data
1857   // Input parameters:
1858   //    deltaT   - smoothing window (in seconds)
1859   //    cutAbs   - max distance of the time info to the median (in time bins)
1860   //    cutSigma - max distance (in the RMS)
1861   //    pcstream - optional debug streamer to store original and filtered info
1862   // Hardwired parameters:
1863   //    kMinPoints =10;       // minimal number of points to define the CE
1864   //    kMinSectors=12;       // minimal number of sectors to define sideCE
1865   // Algorithm:
1866   // 0. Filter almost emty graphs (kMinPoints=10)
1867   // 1. calculate median and RMS per side
1868   // 2. Filter graphs - in respect with side medians
1869   //                  - cutAbs and cutDelta used
1870   // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
1871   // 4. Calculate mean for A side and C side
1872   //
1873   const Int_t kMinPoints =10;       // minimal number of points to define the CE
1874   const Int_t kMinSectors=12;       // minimal number of sectors to define sideCE
1875   const Int_t kMinTime   =400;     // minimal arrival time of CE
1876   TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
1877   Double_t medianY=0;
1878   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1879   if (!cearray) return;
1880   Double_t tmin=-1;
1881   Double_t tmax=-1;
1882   //
1883   //
1884   AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
1885   AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernPressure");
1886   if ( tempMapCE && cavernPressureCE){
1887     //
1888     Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
1889     FilterSensor(cavernPressureCE,960,1050,10, 5.);
1890     if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
1891     if (isOK)  {
1892       // recalculate P/T correction map for time of the CE
1893       AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
1894       driftCalib->SetName("driftPTCE");
1895       driftCalib->SetTitle("driftPTCE");
1896       cearray->AddLast(driftCalib);
1897     }
1898   }
1899   //
1900   // 0. Filter almost emty graphs
1901   //
1902
1903   for (Int_t i=0; i<72;i++){
1904     TGraph *graph= (TGraph*)arrT->At(i);
1905     if (!graph) continue;
1906     if (graph->GetN()<kMinPoints){
1907       arrT->AddAt(0,i);
1908       delete graph;  // delete empty graph
1909       continue;
1910     }
1911     if (tmin<0) tmin = graph->GetX()[0];
1912     if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
1913     //
1914     if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
1915     if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
1916   }
1917   //
1918   // 1. calculate median and RMS per side
1919   //
1920   TArrayF arrA(100000), arrC(100000);
1921   Int_t nA=0, nC=0;
1922   Double_t medianA=0, medianC=0;
1923   Double_t rmsA=0, rmsC=0;
1924   for (Int_t isec=0; isec<72;isec++){
1925     TGraph *graph= (TGraph*)arrT->At(isec);
1926     if (!graph) continue;
1927     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1928       if (graph->GetY()[ipoint]<kMinTime) continue;
1929       if (nA>=arrA.fN) arrA.Set(nA*2);
1930       if (nC>=arrC.fN) arrC.Set(nC*2);
1931       if (isec%36<18)  arrA[nA++]= graph->GetY()[ipoint];
1932       if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
1933     }
1934   }
1935   if (nA>0){
1936     medianA=TMath::Median(nA,arrA.fArray);
1937     rmsA   =TMath::RMS(nA,arrA.fArray);
1938   }
1939   if (nC>0){
1940     medianC=TMath::Median(nC,arrC.fArray);
1941     rmsC   =TMath::RMS(nC,arrC.fArray);
1942   }
1943   //
1944   // 2. Filter graphs - in respect with side medians
1945   //
1946   TArrayD vecX(100000), vecY(100000);
1947   for (Int_t isec=0; isec<72;isec++){
1948     TGraph *graph= (TGraph*)arrT->At(isec);
1949     if (!graph) continue;
1950     Double_t median = (isec%36<18) ? medianA: medianC;
1951     Double_t rms    = (isec%36<18) ? rmsA:    rmsC;
1952     Int_t naccept=0;
1953     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1954       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
1955       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
1956       vecX[naccept]= graph->GetX()[ipoint];
1957       vecY[naccept]= graph->GetY()[ipoint];
1958       naccept++;
1959     }
1960     if (naccept<kMinPoints){
1961       arrT->AddAt(0,isec);
1962       delete graph;  // delete empty graph
1963       continue;
1964     }
1965     TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
1966     delete graph;
1967     arrT->AddAt(graph2,isec);
1968   }
1969   //
1970   // 3. Cut in respect wit the graph median
1971   //
1972   for (Int_t i=0; i<72;i++){
1973     TGraph *graph= (TGraph*)arrT->At(i);
1974     if (!graph) continue;
1975     //
1976     // filter in range
1977     //
1978     TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
1979     if (!graphTS0) continue;
1980     if (graphTS0->GetN()<kMinPoints) {
1981       delete graphTS0;
1982       delete graph;
1983       arrT->AddAt(0,i);
1984       continue;
1985     }
1986     TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
1987     graphTS->Sort();
1988     AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
1989     if (pcstream){
1990       Int_t run = AliTPCcalibDB::Instance()->GetRun();
1991       (*pcstream)<<"filterCE"<<
1992         "run="<<run<<
1993         "isec="<<i<<
1994         "mY="<<medianY<<
1995         "graph.="<<graph<<
1996         "graphTS0.="<<graphTS0<<
1997         "graphTS.="<<graphTS<<
1998         "\n";
1999     }
2000     delete graphTS0;
2001     if (!graphTS) continue;
2002     arrT->AddAt(graphTS,i);
2003     delete graph;
2004   }
2005   //
2006   // Recalculate the mean time A side C side
2007   //
2008   TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2009   Int_t meanPoints=(nA+nC)/72;  // mean number of points
2010   for (Int_t itime=0; itime<200; itime++){
2011     nA=0, nC=0;
2012     Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2013     for (Int_t i=0; i<72;i++){
2014       TGraph *graph= (TGraph*)arrT->At(i);
2015       if (!graph) continue;
2016       if (graph->GetN()<(meanPoints/4)) continue;
2017       if ( (i%36)<18 )  arrA[nA++]=graph->Eval(time);
2018       if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2019     }
2020     xA[itime]=time;
2021     xC[itime]=time;
2022     yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2023     yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2024     eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2025     eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2026   }
2027   //
2028   Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2029   Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2030   if (pcstream){
2031     Int_t run = AliTPCcalibDB::Instance()->GetRun();
2032     (*pcstream)<<"filterAC"<<
2033       "run="<<run<<
2034       "nA="<<nA<<
2035       "nC="<<nC<<
2036       "rmsTA="<<rmsTA<<
2037       "rmsTC="<<rmsTC<<
2038       "\n";
2039   }
2040   //
2041   TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2042   TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2043   TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2044   if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2045   TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2046   if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2047   delete grA;
2048   delete grC;
2049   if (nA<kMinSectors) arrT->AddAt(0,72);
2050   else arrT->AddAt(graphTSA,72);
2051   if (nC<kMinSectors) arrT->AddAt(0,73);
2052   else arrT->AddAt(graphTSC,73);
2053 }
2054
2055
2056 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
2057   //
2058   // Filter Drift velocity measurement using the tracks
2059   // 0.  remove outlyers - error based
2060   //     cutSigma
2061   //
2062   //
2063   const Int_t kMinPoints=1;  // minimal number of points to define value
2064   TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2065   Double_t medianY=0;
2066   if (!arrT) return;
2067   for (Int_t i=0; i<arrT->GetEntries();i++){
2068     TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2069     if (!graph) continue;
2070     if (graph->GetN()<kMinPoints){
2071       delete graph;
2072       arrT->AddAt(0,i);
2073       continue;
2074     }
2075     TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2076     if (!graph2) {
2077       delete graph; arrT->AddAt(0,i); continue;
2078     }
2079     if (graph2->GetN()<1) {
2080       delete graph; arrT->AddAt(0,i); continue;
2081     }
2082     graph2->SetName(graph->GetName());
2083     graph2->SetTitle(graph->GetTitle());
2084     arrT->AddAt(graph2,i);
2085     if (pcstream){
2086       (*pcstream)<<"filterTracks"<<
2087         "run="<<run<<
2088         "isec="<<i<<
2089         "mY="<<medianY<<
2090         "graph.="<<graph<<
2091         "graph2.="<<graph2<<
2092         "\n";
2093     }
2094     delete graph;
2095   }
2096 }
2097
2098
2099 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, TTreeSRedirector *pcstream){
2100   //
2101   // Filter Goofie data
2102   // 0.  remove outlyers      - cutSigma around median
2103   // 1.  smooth the graphs    - deltaT - smoothing window
2104   //
2105   for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2106     //
2107     //
2108     AliDCSSensor *sensor = goofieArray->GetSensor(isensor);
2109     Double_t medianY;
2110     if (sensor &&  sensor->GetGraph()){
2111       TGraph * graph = sensor->GetGraph();
2112       if (isensor==3 && graph->GetN()>1){ // drift velocity
2113         TGraph * graphv = FilterGraphMedianAbs(graph,0.2,medianY);
2114         delete graph;
2115         graph=graphv;
2116       }
2117       TGraph * graph2 = FilterGraphMedian(graph,cutSigma,medianY);
2118       if (!graph2) continue;
2119       AliTPCcalibDButil::SmoothGraph(graph2,deltaT);
2120       if (pcstream){
2121         Int_t run = AliTPCcalibDB::Instance()->GetRun();
2122         (*pcstream)<<"filterG"<<
2123           "run="<<run<<
2124           "isensor="<<isensor<<
2125           "mY="<<medianY<<
2126           "graph.="<<graph<<
2127           "graph2.="<<graph2<<
2128           "\n";
2129       }
2130       if (graph2->GetN()<2) {delete graph2; continue;}
2131       delete graph;
2132       sensor->SetGraph(graph2);
2133     }
2134   }
2135 }
2136
2137
2138
2139
2140
2141
2142
2143 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2144   //
2145   //
2146   // get laser time offset
2147   // median around timeStamp+-deltaT
2148   // QA - chi2 needed for later usage - to be added
2149   //    - currently cut on error
2150   //
2151   Int_t kMinPoints=1;
2152   Double_t kMinDelay=0.01;
2153   Double_t kMinDelayErr=0.0001;
2154
2155   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2156   if (!array) return 0;
2157   TGraphErrors *tlaser=0;
2158   if (array){
2159     if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2160     if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2161   }
2162   if (!tlaser) return 0;
2163   Int_t npoints0= tlaser->GetN();
2164   if (npoints0==0) return 0;
2165   Double_t *xlaser = new Double_t[npoints0];
2166   Double_t *ylaser = new Double_t[npoints0];
2167   Int_t npoints=0;
2168   for (Int_t i=0;i<npoints0;i++){
2169     //printf("%d\n",i);
2170     if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2171     if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2172     xlaser[npoints]=tlaser->GetX()[npoints];
2173     ylaser[npoints]=tlaser->GetY()[npoints];
2174     npoints++;
2175   }
2176   //
2177   //
2178   Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2179   Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2180   //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2181   if (index0<0) index0=0;
2182   if (index1>=npoints-1) index1=npoints-1;
2183   if (index1-index0<kMinPoints) return 0;
2184   //
2185   //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2186     Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2187   delete [] ylaser;
2188   delete [] xlaser;
2189   return mean;
2190 }