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