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