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