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