TPC/AliTPCcalibDButil.cxx

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