TPC/AliTPCcalibDButil.cxx

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