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