TPC/AliTPCcalibDButil.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // Class providing the calculation of derived quantities (mean,rms,fits,...) //
  20 //       of calibration entries                                              //
  21 /*
  22
  23
  24 */
  25 ////////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TMath.h>
  28 #include <TVectorT.h>
  29 #include <TObjArray.h>
  30 #include <TGraph.h>
  31 #include <TFile.h>
  32 #include <TDirectory.h>
  33 #include <TMap.h>
  34 #include <TGraphErrors.h>
  35 #include <AliCDBStorage.h>
  36 #include <AliDCSSensorArray.h>
  37 #include <AliTPCSensorTempArray.h>
  38 #include <AliDCSSensor.h>
  39 #include <AliLog.h>
  40 #include <AliCDBEntry.h>
  41 #include <AliCDBManager.h>
  42 #include <AliCDBId.h>
  43 #include "AliTPCcalibDB.h"
  44 #include "AliTPCCalPad.h"
  45 #include "AliTPCCalROC.h"
  46 #include "AliTPCROC.h"
  47 #include "AliTPCmapper.h"
  48 #include "AliTPCParam.h"
  49 #include "AliTPCCalibRaw.h"
  50 #include "AliTPCPreprocessorOnline.h"
  51 #include "AliTPCdataQA.h"
  52 #include "AliLog.h"
  53 #include "AliTPCcalibDButil.h"
  54 #include "AliTPCCalibVdrift.h"
  55 #include "AliMathBase.h"
  56 #include "AliRelAlignerKalman.h"
  57
  58 const Float_t kAlmost0=1.e-30;
  59
  60 ClassImp(AliTPCcalibDButil)
  61 AliTPCcalibDButil::AliTPCcalibDButil() :
  62   TObject(),
  63   fCalibDB(0),
  64   fPadNoise(0x0),
  65   fPedestals(0x0),
  66   fPulserTmean(0x0),
  67   fPulserTrms(0x0),
  68   fPulserQmean(0x0),
  69   fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
  70   fCETmean(0x0),
  71   fCETrms(0x0),
  72   fCEQmean(0x0),
  73   fALTROMasked(0x0),
  74   fCalibRaw(0x0),
  75   fDataQA(0x0),
  76   fRefMap(0x0),
  77   fCurrentRefMap(0x0),
  78   fRefValidity(""),
  79   fRefPadNoise(0x0),
  80   fRefPedestals(0x0),
  81   fRefPedestalMasked(0x0),
  82   fRefPulserTmean(0x0),
  83   fRefPulserTrms(0x0),
  84   fRefPulserQmean(0x0),
  85   fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
  86   fRefPulserMasked(0x0),
  87   fRefCETmean(0x0),
  88   fRefCETrms(0x0),
  89   fRefCEQmean(0x0),
  90   fRefCEMasked(0x0),
  91   fRefALTROFPED(0x0),
  92   fRefALTROZsThr(0x0),
  93   fRefALTROAcqStart(0x0),
  94   fRefALTROAcqStop(0x0),
  95   fRefALTROMasked(0x0),
  96   fRefCalibRaw(0x0),
  97   fRefDataQA(0x0),
  98   fGoofieArray(0x0),
  99   fMapper(new AliTPCmapper(0x0)),
 100   fNpulserOutliers(-1),
 101   fIrocTimeOffset(0),
 102   fCETmaxLimitAbs(1.5),
 103   fPulTmaxLimitAbs(1.5),
 104   fPulQmaxLimitAbs(5),
 105   fPulQminLimit(11),
 106   fRuns(0),                         // run list with OCDB info
 107   fRunsStart(0),                    // start time for given run
 108   fRunsStop(0)                     // stop time for given run
 109 {
 110   //
 111   // Default ctor
 112   //
 113 }
 114 //_____________________________________________________________________________________
 115 AliTPCcalibDButil::~AliTPCcalibDButil()
 116 {
 117   //
 118   // dtor
 119   //
 120   delete fPulserOutlier;
 121   delete fRefPulserOutlier;
 122   delete fMapper;
 123   if (fRefPadNoise) delete fRefPadNoise;
 124   if (fRefPedestals) delete fRefPedestals;
 125   if (fRefPedestalMasked) delete fRefPedestalMasked;
 126   if (fRefPulserTmean) delete fRefPulserTmean;
 127   if (fRefPulserTrms) delete fRefPulserTrms;
 128   if (fRefPulserQmean) delete fRefPulserQmean;
 129   if (fRefPulserMasked) delete fRefPulserMasked;
 130   if (fRefCETmean) delete fRefCETmean;
 131   if (fRefCETrms) delete fRefCETrms;
 132   if (fRefCEQmean) delete fRefCEQmean;
 133   if (fRefCEMasked) delete fRefCEMasked;
 134   if (fRefALTROFPED) delete fRefALTROFPED;
 135   if (fRefALTROZsThr) delete fRefALTROZsThr;
 136   if (fRefALTROAcqStart) delete fRefALTROAcqStart;
 137   if (fRefALTROAcqStop) delete fRefALTROAcqStop;
 138   if (fRefALTROMasked) delete fRefALTROMasked;
 139   if (fRefCalibRaw) delete fRefCalibRaw;
 140   if (fCurrentRefMap) delete fCurrentRefMap;
 141 }
 142 //_____________________________________________________________________________________
 143 void AliTPCcalibDButil::UpdateFromCalibDB()
 144 {
 145   //
 146   // Update pointers from calibDB
 147   //
 148   if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
 149   fCalibDB->UpdateNonRec();  // load all infromation now
 150   fPadNoise=fCalibDB->GetPadNoise();
 151   fPedestals=fCalibDB->GetPedestals();
 152   fPulserTmean=fCalibDB->GetPulserTmean();
 153   fPulserTrms=fCalibDB->GetPulserTrms();
 154   fPulserQmean=fCalibDB->GetPulserQmean();
 155   fCETmean=fCalibDB->GetCETmean();
 156   fCETrms=fCalibDB->GetCETrms();
 157   fCEQmean=fCalibDB->GetCEQmean();
 158   fALTROMasked=fCalibDB->GetALTROMasked();
 159   fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
 160   fCalibRaw=fCalibDB->GetCalibRaw();
 161   fDataQA=fCalibDB->GetDataQA();
 162   UpdatePulserOutlierMap();
 163 //   SetReferenceRun();
 164 //   UpdateRefDataFromOCDB();
 165 }
 166 //_____________________________________________________________________________________
 167 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
 168                                       Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad * const outCE)
 169 {
 170   //
 171   // Process the CE data for this run
 172   // the return TVectorD arrays contian the results of the fit
 173   // noutliersCE contains the number of pads marked as outliers,
 174   //   not including masked and edge pads
 175   //
 176
 177   //retrieve CE and ALTRO data
 178   if (!fCETmean){
 179     TString fitString(fitFormula);
 180     fitString.ReplaceAll("++","#");
 181     Int_t ndim=fitString.CountChar('#')+2;
 182     fitResultsA.ResizeTo(ndim);
 183     fitResultsC.ResizeTo(ndim);
 184     fitResultsA.Zero();
 185     fitResultsC.Zero();
 186     noutliersCE=-1;
 187     return;
 188   }
 189   noutliersCE=0;
 190   //create outlier map
 191   AliTPCCalPad *out=0;
 192   if (outCE) out=outCE;
 193   else out=new AliTPCCalPad("outCE","outCE");
 194   AliTPCCalROC *rocMasked=0x0;
 195   //loop over all channels
 196   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 197     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 198     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
 199     AliTPCCalROC *rocOut=out->GetCalROC(iroc);
 200     if (!rocData) {
 201       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
 202       rocOut->Add(1.);
 203       continue;
 204     }
 205     //add time offset to IROCs
 206     if (iroc<AliTPCROC::Instance()->GetNInnerSector())
 207       rocData->Add(fIrocTimeOffset);
 208     //select outliers
 209     UInt_t nrows=rocData->GetNrows();
 210     for (UInt_t irow=0;irow<nrows;++irow){
 211       UInt_t npads=rocData->GetNPads(irow);
 212       for (UInt_t ipad=0;ipad<npads;++ipad){
 213         rocOut->SetValue(irow,ipad,0);
 214         //exclude masked pads
 215         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
 216           rocOut->SetValue(irow,ipad,1);
 217           continue;
 218         }
 219         //exclude first two rows in IROC and last two rows in OROC
 220         if (iroc<36){
 221           if (irow<2) rocOut->SetValue(irow,ipad,1);
 222         } else {
 223           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
 224         }
 225         //exclude edge pads
 226         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
 227         Float_t valTmean=rocData->GetValue(irow,ipad);
 228         //exclude values that are exactly 0
 229         if ( !(TMath::Abs(valTmean)>kAlmost0) ) {
 230           rocOut->SetValue(irow,ipad,1);
 231           ++noutliersCE;
 232         }
 233         // exclude channels with too large variations
 234         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
 235           rocOut->SetValue(irow,ipad,1);
 236           ++noutliersCE;
 237         }
 238       }
 239     }
 240   }
 241   //perform fit
 242   TMatrixD dummy;
 243   Float_t chi2Af,chi2Cf;
 244   fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
 245   chi2A=chi2Af;
 246   chi2C=chi2Cf;
 247   if (!outCE) delete out;
 248 }
 249 //_____________________________________________________________________________________
 250 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
 251                      TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
 252                      Float_t &driftTimeA, Float_t &driftTimeC )
 253 {
 254   //
 255   // Calculate statistical information from the CE graphs for drift time and charge
 256   //
 257
 258   //reset arrays
 259   vecTEntries.ResizeTo(72);
 260   vecTMean.ResizeTo(72);
 261   vecTRMS.ResizeTo(72);
 262   vecTMedian.ResizeTo(72);
 263   vecQEntries.ResizeTo(72);
 264   vecQMean.ResizeTo(72);
 265   vecQRMS.ResizeTo(72);
 266   vecQMedian.ResizeTo(72);
 267   vecTEntries.Zero();
 268   vecTMean.Zero();
 269   vecTRMS.Zero();
 270   vecTMedian.Zero();
 271   vecQEntries.Zero();
 272   vecQMean.Zero();
 273   vecQRMS.Zero();
 274   vecQMedian.Zero();
 275   driftTimeA=0;
 276   driftTimeC=0;
 277   TObjArray *arrT=fCalibDB->GetCErocTtime();
 278   TObjArray *arrQ=fCalibDB->GetCErocQtime();
 279   if (arrT){
 280     for (Int_t isec=0;isec<74;++isec){
 281       TGraph *gr=(TGraph*)arrT->At(isec);
 282       if (!gr) continue;
 283       TVectorD values;
 284       Int_t npoints = gr->GetN();
 285       values.ResizeTo(npoints);
 286       Int_t nused =0;
 287       //skip first points, theres always some problems with finding the CE position
 288       for (Int_t ipoint=4; ipoint<npoints; ipoint++){
 289         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
 290           values[nused]=gr->GetY()[ipoint];
 291           nused++;
 292         }
 293       }
 294       //
 295       if (isec<72) vecTEntries[isec]= nused;
 296       if (nused>1){
 297         if (isec<72){
 298           vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 299           vecTMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 300           vecTRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 301         } else if (isec==72){
 302           driftTimeA=TMath::Median(nused,values.GetMatrixArray());
 303         } else if (isec==73){
 304           driftTimeC=TMath::Median(nused,values.GetMatrixArray());
 305         }
 306       }
 307     }
 308   }
 309   if (arrQ){
 310     for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
 311       TGraph *gr=(TGraph*)arrQ->At(isec);
 312       if (!gr) continue;
 313       TVectorD values;
 314       Int_t npoints = gr->GetN();
 315       values.ResizeTo(npoints);
 316       Int_t nused =0;
 317       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 318         if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
 319           values[nused]=gr->GetY()[ipoint];
 320           nused++;
 321         }
 322       }
 323       //
 324       vecQEntries[isec]= nused;
 325       if (nused>1){
 326         vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 327         vecQMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 328         vecQRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 329       }
 330     }
 331   }
 332 }
 333
 334 //_____________________________________________________________________________________
 335 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
 336                       TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
 337                       Int_t &nonMaskedZero, Int_t &nNaN)
 338 {
 339   //
 340   // process noise data
 341   // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
 342   //    OROCs small pads [2] and OROCs large pads [3]
 343   // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
 344   // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
 345   //
 346
 347   //set proper size and reset
 348   const UInt_t infoSize=4;
 349   vNoiseMean.ResizeTo(infoSize);
 350   vNoiseMeanSenRegions.ResizeTo(infoSize);
 351   vNoiseRMS.ResizeTo(infoSize);
 352   vNoiseRMSSenRegions.ResizeTo(infoSize);
 353   vNoiseMean.Zero();
 354   vNoiseMeanSenRegions.Zero();
 355   vNoiseRMS.Zero();
 356   vNoiseRMSSenRegions.Zero();
 357   nonMaskedZero=0;
 358   nNaN=0;
 359   //counters
 360   TVectorD c(infoSize);
 361   TVectorD cs(infoSize);
 362   //tpc parameters
 363   AliTPCParam par;
 364   par.Update();
 365   //retrieve noise and ALTRO data
 366   if (!fPadNoise) return;
 367   AliTPCCalROC *rocMasked=0x0;
 368   //create IROC, OROC1, OROC2 and sensitive region masks
 369   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 370     AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
 371     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 372     UInt_t nrows=noiseROC->GetNrows();
 373     for (UInt_t irow=0;irow<nrows;++irow){
 374       UInt_t npads=noiseROC->GetNPads(irow);
 375       for (UInt_t ipad=0;ipad<npads;++ipad){
 376         //don't use masked channels;
 377         if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
 378         Float_t noiseVal=noiseROC->GetValue(irow,ipad);
 379         //check if noise==0
 380         if (noiseVal<kAlmost0) {
 381           ++nonMaskedZero;
 382           continue;
 383         }
 384         //check for nan
 385         if ( !(noiseVal<10000000) ){
 386 //           printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
 387           ++nNaN;
 388           continue;
 389         }
 390         Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
 391         Int_t masksen=1; // sensitive pards are not masked (0)
 392         if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
 393         if (isec<AliTPCROC::Instance()->GetNInnerSector()){
 394           //IROCs
 395           if (irow>19&&irow<46){
 396             if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
 397           }
 398           Int_t type=1;
 399           vNoiseMean[type]+=noiseVal;
 400           vNoiseRMS[type]+=noiseVal*noiseVal;
 401           ++c[type];
 402           if (!masksen){
 403             vNoiseMeanSenRegions[type]+=noiseVal;
 404             vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 405             ++cs[type];
 406           }
 407         } else {
 408           //OROCs
 409           //define sensive regions
 410           if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
 411           if ( irow>75 ){
 412             Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
 413             if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
 414           }
 415           if ((Int_t)irow<par.GetNRowUp1()){
 416             //OROC1
 417             Int_t type=2;
 418             vNoiseMean[type]+=noiseVal;
 419             vNoiseRMS[type]+=noiseVal*noiseVal;
 420             ++c[type];
 421             if (!masksen){
 422               vNoiseMeanSenRegions[type]+=noiseVal;
 423               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 424               ++cs[type];
 425             }
 426           }else{
 427             //OROC2
 428             Int_t type=3;
 429             vNoiseMean[type]+=noiseVal;
 430             vNoiseRMS[type]+=noiseVal*noiseVal;
 431             ++c[type];
 432             if (!masksen){
 433               vNoiseMeanSenRegions[type]+=noiseVal;
 434               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 435               ++cs[type];
 436             }
 437           }
 438         }
 439         //whole tpc
 440         Int_t type=0;
 441         vNoiseMean[type]+=noiseVal;
 442         vNoiseRMS[type]+=noiseVal*noiseVal;
 443         ++c[type];
 444         if (!masksen){
 445           vNoiseMeanSenRegions[type]+=noiseVal;
 446           vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 447           ++cs[type];
 448         }
 449       }//end loop pads
 450     }//end loop rows
 451   }//end loop sectors (rocs)
 452
 453   //calculate mean and RMS
 454   const Double_t verySmall=0.0000000001;
 455   for (UInt_t i=0;i<infoSize;++i){
 456     Double_t mean=0;
 457     Double_t rms=0;
 458     Double_t meanSen=0;
 459     Double_t rmsSen=0;
 460
 461     if (c[i]>verySmall){
 462 //       printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
 463       mean=vNoiseMean[i]/c[i];
 464       rms=vNoiseRMS[i];
 465       rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
 466     }
 467     vNoiseMean[i]=mean;
 468     vNoiseRMS[i]=rms;
 469
 470     if (cs[i]>verySmall){
 471       meanSen=vNoiseMeanSenRegions[i]/cs[i];
 472       rmsSen=vNoiseRMSSenRegions[i];
 473       rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
 474     }
 475     vNoiseMeanSenRegions[i]=meanSen;
 476     vNoiseRMSSenRegions[i]=rmsSen;
 477   }
 478 }
 479
 480 //_____________________________________________________________________________________
 481 void AliTPCcalibDButil::ProcessQAData(TVectorD &vQaOcc, TVectorD &vQaQtot,
 482                                       TVectorD &vQaQmax)
 483 {
 484   //
 485   // process QA data
 486   //
 487   // vQaOcc/Qtot/Qmax contains the Mean occupancy/Qtot/Qmax for each sector
 488   //
 489
 490
 491   const UInt_t infoSize = 72;
 492   //reset counters to error number
 493   vQaOcc.ResizeTo(infoSize);
 494   vQaOcc.Zero();
 495   vQaQtot.ResizeTo(infoSize);
 496   vQaQtot.Zero();
 497   vQaQmax.ResizeTo(infoSize);
 498   vQaQmax.Zero();
 499   //counter
 500   //retrieve pulser and ALTRO data
 501
 502   if (!fDataQA) {
 503
 504     AliInfo("No QA data");
 505     return;
 506   }
 507   if (fDataQA->GetEventCounter()<=0) {
 508
 509     AliInfo("No QA data");
 510     return; // no data processed
 511   }
 512   //
 513   fDataQA->Analyse();
 514
 515   TVectorD normOcc(infoSize);
 516   TVectorD normQ(infoSize);
 517
 518   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 519
 520     printf("Sector %d\n", isec);
 521
 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     printf("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(pq/pqRef-1);
 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 (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
1405     AliTPCCalROC *rocOut       = out->GetCalROC(iroc);
1406     AliTPCCalROC *rocCEQ       = fCEQmean->GetCalROC(iroc);
1407     AliTPCCalROC *rocCETrms    = fCETrms->GetCalROC(iroc);
1408     Double_t trocMedian        = rocData->GetMedian();
1409     //
1410     if (!rocData) {
1411       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1412       rocOut->Add(1.);
1413       continue;
1414     }
1415     //
1416     //select outliers
1417     UInt_t nrows=rocData->GetNrows();
1418     for (UInt_t irow=0;irow<nrows;++irow){
1419       UInt_t npads=rocData->GetNPads(irow);
1420       for (UInt_t ipad=0;ipad<npads;++ipad){
1421         rocOut->SetValue(irow,ipad,0);
1422         Float_t valTmean=rocData->GetValue(irow,ipad);
1423         Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1424         Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1425         //0. exclude masked pads
1426         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1427           rocOut->SetValue(irow,ipad,1);
1428           continue;
1429         }
1430         //1. exclude first two rows in IROC and last two rows in OROC
1431         if (iroc<36){
1432           if (irow<2) rocOut->SetValue(irow,ipad,1);
1433         } else {
1434           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1435         }
1436         //2. exclude edge pads
1437         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1438         //exclude values that are exactly 0
1439         if ( TMath::Abs(valTmean)<kAlmost0) {
1440           rocOut->SetValue(irow,ipad,1);
1441           ++noutliersCE;
1442         }
1443         //3.  exclude channels with too large variations
1444         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1445           rocOut->SetValue(irow,ipad,1);
1446           ++noutliersCE;
1447         }
1448         //
1449         //4.  exclude channels with too small signal
1450         if (valQmean<minSignal) {
1451           rocOut->SetValue(irow,ipad,1);
1452           ++noutliersCE;
1453         }
1454         //
1455         //5. exclude channels with too small rms
1456         if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1457           rocOut->SetValue(irow,ipad,1);
1458           ++noutliersCE;
1459         }
1460         //
1461         //6. exclude channels to far from the chamber median
1462         if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1463           rocOut->SetValue(irow,ipad,1);
1464           ++noutliersCE;
1465         }
1466       }
1467     }
1468   }
1469   //
1470   return out;
1471 }
1472
1473
1474 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad * const pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1475   //
1476   // Author: marian.ivanov@cern.ch
1477   //
1478   // Create outlier map for Pulser
1479   // Parameters:
1480   //  Return value     - outlyer map
1481   //  noutlyersPulser  - number of outlyers
1482   //  cutTime          - absolute cut - distance to the median of chamber
1483   //  cutnRMSQ         - nsigma cut from median  q distribution per chamber
1484   //  cutnRMSrms       - nsigma cut from median  rms distribution
1485   // Outlyers criteria:
1486   // 0. Exclude masked pads
1487   // 1. Exclude time outlyers (default 3 time bins)
1488   // 2. Exclude q outlyers    (default 5 sigma)
1489   // 3. Exclude rms outlyers  (default 5 sigma)
1490   noutliersPulser=0;
1491   AliTPCCalPad *out=pulserOut;
1492   if (!out)     out= new AliTPCCalPad("outPulser","outPulser");
1493   AliTPCCalROC *rocMasked=0x0;
1494   if (!fPulserTmean) return 0;
1495   if (!fPulserTrms) return 0;
1496   if (!fPulserQmean) return 0;
1497   //
1498   //loop over all channels
1499   //
1500   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1501     if (fALTROMasked)   rocMasked= fALTROMasked->GetCalROC(iroc);
1502     AliTPCCalROC *rocData       = fPulserTmean->GetCalROC(iroc);
1503     AliTPCCalROC *rocOut        = out->GetCalROC(iroc);
1504     AliTPCCalROC *rocPulserQ    = fPulserQmean->GetCalROC(iroc);
1505     AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1506     //
1507     Double_t rocMedianT         = rocData->GetMedian();
1508     Double_t rocMedianQ         = rocPulserQ->GetMedian();
1509     Double_t rocRMSQ            = rocPulserQ->GetRMS();
1510     Double_t rocMedianTrms      = rocPulserTrms->GetMedian();
1511     Double_t rocRMSTrms         = rocPulserTrms->GetRMS();
1512     for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1513       rocOut->SetValue(ichannel,0);
1514       Float_t valTmean=rocData->GetValue(ichannel);
1515       Float_t valQmean=rocPulserQ->GetValue(ichannel);
1516       Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1517       Int_t isOut=0;
1518       if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1519       if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1520       if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1521       rocOut->SetValue(ichannel,isOut);
1522       if (isOut) noutliersPulser++;
1523     }
1524   }
1525   return out;
1526 }
1527
1528
1529 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1530   //
1531   // Author : Marian Ivanov
1532   // Create pad time0 correction map using information from the CE and from pulser
1533   //
1534   //
1535   // Return PadTime0 to be used for time0 relative alignment
1536   // if dump file specified intermediat results are dumped to the fiel and can be visualized
1537   // using $ALICE_ROOT/TPC/script/gui application
1538   //
1539   // fitResultsA - fitParameters A side
1540   // fitResultsC - fitParameters C side
1541   // chi2A       - chi2/ndf for A side (assuming error 1 time bin)
1542   // chi2C       - chi2/ndf for C side (assuming error 1 time bin)
1543   //
1544   //
1545   // Algorithm:
1546   // 1. Find outlier map for CE
1547   // 2. Find outlier map for Pulser
1548   // 3. Replace outlier by median at given sector  (median without outliers)
1549   // 4. Substract from the CE data pulser
1550   // 5. Fit the CE with formula
1551   //    5.1) (IROC-OROC) offset
1552   //    5.2) gx
1553   //    5.3) gy
1554   //    5.4) (lx-xmid)
1555   //    5.5) (IROC-OROC)*(lx-xmid)
1556   //    5.6) (ly/lx)^2
1557   // 6. Substract gy fit dependence from the CE data
1558   // 7. Add pulser back to CE data
1559   // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1560   // 9. return CE data
1561   //
1562   // Time0 <= padCE = padCEin  -padCEfitGy  - if not outlier
1563   // Time0 <= padCE = padFitAll-padCEfitGy  - if outlier
1564
1565   // fit formula
1566   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)";
1567   // output for fit formula
1568   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)";
1569   // gy part of formula
1570   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)";
1571   //
1572   //
1573   if (!fCETmean) return 0;
1574   Double_t pgya,pgyc,pchi2a,pchi2c;
1575   AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1576   AliTPCCalPad * padCEOut     = CreateCEOutlyerMap(nOut);
1577
1578   AliTPCCalPad * padPulser    = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1579   AliTPCCalPad * padCE        = new AliTPCCalPad(*fCETmean);
1580   AliTPCCalPad * padCEIn      = new AliTPCCalPad(*fCETmean);
1581   AliTPCCalPad * padOut       = new AliTPCCalPad("padOut","padOut");
1582   padPulser->SetName("padPulser");
1583   padPulserOut->SetName("padPulserOut");
1584   padCE->SetName("padCE");
1585   padCEIn->SetName("padCEIn");
1586   padCEOut->SetName("padCEOut");
1587   padOut->SetName("padOut");
1588
1589   //
1590   // make combined outlyers map
1591   // and replace outlyers in maps with median for chamber
1592   //
1593   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1594     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1595     AliTPCCalROC * rocPulser    = padPulser->GetCalROC(iroc);
1596     AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1597     AliTPCCalROC * rocCEOut     = padCEOut->GetCalROC(iroc);
1598     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1599     Double_t ceMedian           = rocCE->GetMedian(rocCEOut);
1600     Double_t pulserMedian       = rocPulser->GetMedian(rocCEOut);
1601     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1602       if (rocPulserOut->GetValue(ichannel)>0) {
1603         rocPulser->SetValue(ichannel,pulserMedian);
1604         rocOut->SetValue(ichannel,1);
1605       }
1606       if (rocCEOut->GetValue(ichannel)>0) {
1607         rocCE->SetValue(ichannel,ceMedian);
1608         rocOut->SetValue(ichannel,1);
1609       }
1610     }
1611   }
1612   //
1613   // remove pulser time 0
1614   //
1615   padCE->Add(padPulser,-1);
1616   //
1617   // Make fits
1618   //
1619   TMatrixD dummy;
1620   Float_t chi2Af,chi2Cf;
1621   padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1622   chi2A=chi2Af;
1623   chi2C=chi2Cf;
1624   //
1625   AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1626   padCEFitGY->SetName("padCEFitGy");
1627   //
1628   AliTPCCalPad *padCEFit  =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1629   padCEFit->SetName("padCEFit");
1630   //
1631   AliTPCCalPad* padCEDiff  = new AliTPCCalPad(*padCE);
1632   padCEDiff->SetName("padCEDiff");
1633   padCEDiff->Add(padCEFit,-1.);
1634   //
1635   //
1636   padCE->Add(padCEFitGY,-1.);
1637
1638   padCE->Add(padPulser,1.);
1639   Double_t padmedian = padCE->GetMedian();
1640   padCE->Add(-padmedian);  // normalize to median
1641   //
1642   // Replace outliers by fit value - median of diff per given chamber -GY fit
1643   //
1644   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1645     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1646     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1647     AliTPCCalROC * rocCEFit     = padCEFit->GetCalROC(iroc);
1648     AliTPCCalROC * rocCEFitGY   = padCEFitGY->GetCalROC(iroc);
1649     AliTPCCalROC * rocCEDiff    = padCEDiff->GetCalROC(iroc);
1650     //
1651     Double_t diffMedian         = rocCEDiff->GetMedian(rocOut);
1652     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1653       if (rocOut->GetValue(ichannel)==0) continue;
1654       Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1655       rocCE->SetValue(ichannel,value);
1656     }
1657   }
1658   //
1659   //
1660   if (dumpfile){
1661     //dump to the file - result can be visualized
1662     AliTPCPreprocessorOnline preprocesor;
1663     preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1664     preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1665     preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1666     preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1667     //
1668     preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1669     preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1670     //
1671     preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1672     preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1673     preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1674     preprocesor.DumpToFile(dumpfile);
1675   }
1676   delete padPulser;
1677   delete padPulserOut;
1678   delete padCEIn;
1679   delete padCEOut;
1680   delete padOut;
1681   delete padCEDiff;
1682   delete padCEFitGY;
1683   return padCE;
1684 }
1685
1686
1687
1688
1689
1690 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1691   //
1692   // find the closest point to xref  in x  direction
1693   // return dx and value
1694   dx = 0;
1695   y = 0;
1696
1697   if(!graph) return 0;
1698   if(graph->GetN() < 1) return 0;
1699
1700   Int_t index=0;
1701   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1702   if (index<0) index=0;
1703   if(graph->GetN()==1) {
1704     dx = xref-graph->GetX()[index];
1705   }
1706   else {
1707     if (index>=graph->GetN()-1) index=graph->GetN()-2;
1708     if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1709     dx = xref-graph->GetX()[index];
1710   }
1711   y  = graph->GetY()[index];
1712   return index;
1713 }
1714
1715 Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1716   //
1717   // Get the correction of the trigger offset
1718   // combining information from the laser track calibration
1719   // and from cosmic calibration
1720   //
1721   // run       - run number
1722   // timeStamp - tim stamp in seconds
1723   // deltaT    - integration period to calculate offset
1724   // deltaTLaser -max validity of laser data
1725   // valType   - 0 - median, 1- mean
1726   //
1727   // Integration vaues are just recomendation - if not possible to get points
1728   // automatically increase the validity by factor 2
1729   // (recursive algorithm until one month of data taking)
1730   //
1731   //
1732   const Float_t kLaserCut=0.0005;
1733   const Int_t   kMaxPeriod=3600*24*30*12; // one year max
1734   const Int_t   kMinPoints=20;
1735   //
1736   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1737   if (!array) {
1738     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1739   }
1740   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1741   if (!array) return 0;
1742   //
1743   TGraphErrors *laserA[3]={0,0,0};
1744   TGraphErrors *laserC[3]={0,0,0};
1745   TGraphErrors *cosmicAll=0;
1746   laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1747   laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1748   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1749   //
1750   //
1751   if (!cosmicAll) return 0;
1752   Int_t nmeasC=cosmicAll->GetN();
1753   Float_t *tdelta = new Float_t[nmeasC];
1754   Int_t nused=0;
1755   for (Int_t i=0;i<nmeasC;i++){
1756     if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1757     Float_t ccosmic=cosmicAll->GetY()[i];
1758     Double_t yA=0,yC=0,dA=0,dC=0;
1759     if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1760     if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1761     //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1762     //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1763     //
1764     if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1765     Float_t claser=0;
1766     if (TMath::Abs(yA-yC)<kLaserCut) {
1767       claser=(yA-yC)*0.5;
1768     }else{
1769       if (i%2==0)  claser=yA;
1770       if (i%2==1)  claser=yC;
1771     }
1772     tdelta[nused]=ccosmic-claser;
1773     nused++;
1774   }
1775   if (nused<kMinPoints &&deltaT<kMaxPeriod) return  AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1776   if (nused<kMinPoints) {
1777     printf("AliFatal: No time offset calibration available\n");
1778     return 0;
1779   }
1780   Double_t median = TMath::Median(nused,tdelta);
1781   Double_t mean  = TMath::Mean(nused,tdelta);
1782   delete [] tdelta;
1783   return (valType==0) ? median:mean;
1784 }
1785
1786 Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1787   //
1788   // Get the correction of the drift velocity
1789   // combining information from the laser track calibration
1790   // and from cosmic calibration
1791   //
1792   // dist      - return value - distance to closest point in graph
1793   // run       - run number
1794   // timeStamp - tim stamp in seconds
1795   // deltaT    - integration period to calculate time0 offset
1796   // deltaTLaser -max validity of laser data
1797   // valType   - 0 - median, 1- mean
1798   //
1799   // Integration vaues are just recomendation - if not possible to get points
1800   // automatically increase the validity by factor 2
1801   // (recursive algorithm until one month of data taking)
1802   //
1803   //
1804   //
1805   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1806   if (!array) {
1807     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1808   }
1809   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1810   if (!array) return 0;
1811   TGraphErrors *cosmicAll=0;
1812   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1813   if (!cosmicAll) return 0;
1814   Double_t grY=0;
1815   AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1816
1817   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1818   Double_t vcosmic =  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1819   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1820   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
1821   return  vcosmic-t0;
1822
1823   /*
1824     Example usage:
1825
1826     Int_t run=89000
1827     TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1828     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1829     laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1830     //
1831     Double_t *yvd= new Double_t[cosmicAll->GetN()];
1832     Double_t *yt0= new Double_t[cosmicAll->GetN()];
1833     for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1834     for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1835
1836     TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1837     TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1838
1839   */
1840
1841 }
1842
1843 const char* AliTPCcalibDButil::GetGUIRefTreeDefaultName()
1844 {
1845   //
1846   // Create a default name for the gui file
1847   //
1848
1849   return Form("guiRefTreeRun%s.root",GetRefValidity());
1850 }
1851
1852 Bool_t AliTPCcalibDButil::CreateGUIRefTree(const char* filename)
1853 {
1854   //
1855   // Create a gui reference tree
1856   // if dirname and filename are empty default values will be used
1857   // this is the recommended way of using this function
1858   // it allows to check whether a file with the given run validity alredy exists
1859   //
1860   if (!AliCDBManager::Instance()->GetDefaultStorage()){
1861     AliError("Default Storage not set. Cannot create reference calibration Tree!");
1862     return kFALSE;
1863   }
1864
1865   TString file=filename;
1866   if (file.IsNull()) file=GetGUIRefTreeDefaultName();
1867
1868   AliTPCPreprocessorOnline prep;
1869   //noise and pedestals
1870   if (fRefPedestals) prep.AddComponent(new AliTPCCalPad(*(fRefPedestals)));
1871   if (fRefPadNoise ) prep.AddComponent(new AliTPCCalPad(*(fRefPadNoise)));
1872   if (fRefPedestalMasked) prep.AddComponent(new AliTPCCalPad(*fRefPedestalMasked));
1873   //pulser data
1874   if (fRefPulserTmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTmean)));
1875   if (fRefPulserTrms ) prep.AddComponent(new AliTPCCalPad(*(fRefPulserTrms)));
1876   if (fRefPulserQmean) prep.AddComponent(new AliTPCCalPad(*(fRefPulserQmean)));
1877   if (fRefPulserMasked) prep.AddComponent(new AliTPCCalPad(*fRefPulserMasked));
1878   //CE data
1879   if (fRefCETmean) prep.AddComponent(new AliTPCCalPad(*(fRefCETmean)));
1880   if (fRefCETrms ) prep.AddComponent(new AliTPCCalPad(*(fRefCETrms)));
1881   if (fRefCEQmean) prep.AddComponent(new AliTPCCalPad(*(fRefCEQmean)));
1882   if (fRefCEMasked) prep.AddComponent(new AliTPCCalPad(*fRefCEMasked));
1883   //Altro data
1884   if (fRefALTROAcqStart ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStart )));
1885   if (fRefALTROZsThr    ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROZsThr    )));
1886   if (fRefALTROFPED     ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROFPED     )));
1887   if (fRefALTROAcqStop  ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROAcqStop  )));
1888   if (fRefALTROMasked   ) prep.AddComponent(new AliTPCCalPad(*(fRefALTROMasked   )));
1889   //QA
1890   AliTPCdataQA *dataQA=fRefDataQA;
1891   if (dataQA) {
1892     if (dataQA->GetNLocalMaxima())
1893       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNLocalMaxima())));
1894     if (dataQA->GetMaxCharge())
1895       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMaxCharge())));
1896     if (dataQA->GetMeanCharge())
1897       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetMeanCharge())));
1898     if (dataQA->GetNoThreshold())
1899       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNoThreshold())));
1900     if (dataQA->GetNTimeBins())
1901       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNTimeBins())));
1902     if (dataQA->GetNPads())
1903       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetNPads())));
1904     if (dataQA->GetTimePosition())
1905       prep.AddComponent(new AliTPCCalPad(*(dataQA->GetTimePosition())));
1906   }
1907   prep.DumpToFile(file.Data());
1908   return kTRUE;
1909 }
1910
1911 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1912   //
1913   // Get the correction of the drift velocity using the offline laser tracks calbration
1914   //
1915   // run       - run number
1916   // timeStamp - tim stamp in seconds
1917   // deltaT    - integration period to calculate time0 offset
1918   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1919   // Note in case no data form both A and C side - the value from active side used
1920   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1921
1922   return GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1923 }
1924
1925 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracksOnline(Double_t &dist, Int_t /*run*/, Int_t timeStamp, Double_t deltaT, Int_t side){
1926   //
1927   // Get the correction of the drift velocity using the online laser tracks calbration
1928   //
1929   // run       - run number
1930   // timeStamp - tim stamp in seconds
1931   // deltaT    - integration period to calculate time0 offset
1932   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1933   // Note in case no data form both A and C side - the value from active side used
1934   TObjArray *array =AliTPCcalibDB::Instance()->GetCEfitsDrift();
1935
1936   Double_t dv = GetVDriftTPCLaserTracksCommon(dist, timeStamp, deltaT, side, array);
1937   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
1938   if (!param) return 0;
1939
1940   //the drift velocity is hard wired in the AliTPCCalibCE class, since online there is no access to OCDB
1941   dv*=param->GetDriftV()/2.61301900000000000e+06;
1942   if (dv>1e-20) dv=1/dv-1;
1943   else return 0;
1944   // T/P correction
1945   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1946
1947   AliTPCSensorTempArray *temp = (AliTPCSensorTempArray*)cearray->FindObject("TempMap");
1948   AliDCSSensor *press         = (AliDCSSensor*)cearray->FindObject("CavernAtmosPressure");
1949
1950   Double_t corrPTA=0;
1951   Double_t corrPTC=0;
1952
1953   if (temp&&press) {
1954     AliTPCCalibVdrift corr(temp,press,0);
1955     corrPTA=corr.GetPTRelative(timeStamp,0);
1956     corrPTC=corr.GetPTRelative(timeStamp,1);
1957   }
1958
1959   if (side==0) dv -=  corrPTA;
1960   if (side==1) dv -=  corrPTC;
1961   if (side==2) dv -=  (corrPTA+corrPTC)/2;
1962
1963   return dv;
1964 }
1965
1966 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracksCommon(Double_t &dist, Int_t timeStamp, Double_t deltaT,
1967   Int_t side, TObjArray * const array){
1968   //
1969   // common drift velocity retrieval for online and offline method
1970   //
1971   TGraphErrors *grlaserA=0;
1972   TGraphErrors *grlaserC=0;
1973   Double_t vlaserA=0, vlaserC=0;
1974   if (!array) return 0;
1975   grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1976   grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1977   Double_t deltaY;
1978   if (grlaserA && grlaserA->GetN()>0) {
1979     AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1980     if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
1981     else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1982   }
1983   if (grlaserC && grlaserC->GetN()>0) {
1984     AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1985     if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
1986     else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1987   }
1988   if (side==0) return vlaserA;
1989   if (side==1) return vlaserC;
1990   Double_t mdrift=(vlaserA+vlaserC)*0.5;
1991   if (!grlaserA) return vlaserC;
1992   if (!grlaserC) return vlaserA;
1993   return mdrift;
1994 }
1995
1996
1997 Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1998   //
1999   // Get the correction of the drift velocity using the CE laser data
2000   // combining information from the CE,  laser track calibration
2001   // and P/T calibration
2002   //
2003   // run       - run number
2004   // timeStamp - tim stamp in seconds
2005   // deltaT    - integration period to calculate time0 offset
2006   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
2007   TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
2008   if (!arrT) return 0;
2009   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
2010   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
2011   AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
2012   //
2013   //
2014   Double_t corrPTA = 0, corrPTC=0;
2015   Double_t ltime0A = 0, ltime0C=0;
2016   Double_t gry=0;
2017   Double_t corrA=0, corrC=0;
2018   Double_t timeA=0, timeC=0;
2019   const Double_t kEpsilon = 0.00001;
2020   TGraph *graphA = (TGraph*)arrT->At(72);
2021   TGraph *graphC = (TGraph*)arrT->At(73);
2022   if (!graphA && !graphC) return 0.;
2023   if (graphA &&graphA->GetN()>0) {
2024     AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
2025     timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
2026     Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
2027     ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2028     if(ltime0A < kEpsilon) return 0;
2029     if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
2030     corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2031     corrA-=corrPTA;
2032   }
2033   if (graphC&&graphC->GetN()>0){
2034     AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
2035     timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
2036     Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
2037     ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
2038     if(ltime0C < kEpsilon) return 0;
2039 if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
2040     corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
2041     corrC-=corrPTC;
2042   }
2043
2044   if (side ==0 ) return corrA;
2045   if (side ==1 ) return corrC;
2046   Double_t corrM= (corrA+corrC)*0.5;
2047   if (!graphA) corrM=corrC;
2048   if (!graphC) corrM=corrA;
2049   return corrM;
2050 }
2051
2052 Double_t  AliTPCcalibDButil::GetVDriftTPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2053   //
2054   // return drift velocity using the TPC-ITS matchin method
2055   // return also distance to the closest point
2056   //
2057   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2058   TGraphErrors *graph=0;
2059   dist=0;
2060   if (!array) return 0;
2061   //array->ls();
2062   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSB_TPC_DRIFTVD");
2063   if (!graph) return 0;
2064   Double_t deltaY;
2065   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2066   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2067   return value;
2068 }
2069
2070 Double_t AliTPCcalibDButil::GetTime0TPCITS(Double_t &dist, Int_t run, Int_t timeStamp){
2071   //
2072   // Get time dependent time 0 (trigger delay in cm) correction
2073   // Arguments:
2074   // timestamp - timestamp
2075   // run       - run number
2076   //
2077   // Notice - Extrapolation outside of calibration range  - using constant function
2078   //
2079   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2080   TGraphErrors *graph=0;
2081   dist=0;
2082   if (!array) return 0;
2083   graph = (TGraphErrors*)array->FindObject("ALIGN_ITSM_TPC_T0");
2084   if (!graph) return 0;
2085   Double_t deltaY;
2086   AliTPCcalibDButil::GetNearest(graph,timeStamp,dist,deltaY);
2087   Double_t value = AliTPCcalibDButil::EvalGraphConst(graph,timeStamp);
2088   return value;
2089 }
2090
2091
2092
2093
2094
2095 Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
2096   //
2097   // VERY obscure method - we need something in framework
2098   // Find the TPC runs with temperature OCDB entry
2099   // cache the start and end of the run
2100   //
2101   AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
2102   if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
2103   if (!storage) return 0;
2104   TString path=storage->GetURI();
2105   TString runsT;
2106   {
2107     TString command;
2108     if (path.Contains("local")){  // find the list if local system
2109       path.ReplaceAll("local://","");
2110       path+="TPC/Calib/Temperature";
2111       command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
2112     }
2113     runsT=gSystem->GetFromPipe(command);
2114   }
2115   TObjArray *arr= runsT.Tokenize("r");
2116   if (!arr) return 0;
2117   //
2118   TArrayI indexes(arr->GetEntries());
2119   TArrayI runs(arr->GetEntries());
2120   Int_t naccept=0;
2121   {for (Int_t irun=0;irun<arr->GetEntries();irun++){
2122       Int_t irunN = atoi(arr->At(irun)->GetName());
2123       if (irunN<startRun) continue;
2124       if (irunN>stopRun) continue;
2125       runs[naccept]=irunN;
2126       naccept++;
2127     }}
2128   fRuns.Set(naccept);
2129   fRunsStart.Set(fRuns.fN);
2130   fRunsStop.Set(fRuns.fN);
2131   TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
2132   for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
2133
2134   //
2135   AliCDBEntry * entry = 0;
2136   {for (Int_t irun=0;irun<fRuns.fN; irun++){
2137       entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
2138       if (!entry) continue;
2139       AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
2140       if (!tmpRun) continue;
2141       fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
2142       fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
2143       //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
2144     }}
2145   return fRuns.fN;
2146 }
2147
2148
2149 Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
2150   //
2151   // binary search - find the run for given time stamp
2152   //
2153   Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
2154   Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
2155   Int_t cindex  = -1;
2156   for (Int_t index=index0; index<=index1; index++){
2157     if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
2158     if (debug) {
2159       printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
2160     }
2161   }
2162   if (cindex<0) cindex =(index0+index1)/2;
2163   if (cindex<0) {
2164     return 0;
2165   }
2166   return fRuns[cindex];
2167 }
2168
2169
2170
2171
2172
2173 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
2174   //
2175   // filter outlyer measurement
2176   // Only points around median +- sigmaCut filtered
2177   if (!graph) return  0;
2178   Int_t kMinPoints=2;
2179   Int_t npoints0 = graph->GetN();
2180   Int_t npoints=0;
2181   Float_t  rmsY=0;
2182   //
2183   //
2184   if (npoints0<kMinPoints) return 0;
2185
2186   Double_t *outx=new Double_t[npoints0];
2187   Double_t *outy=new Double_t[npoints0];
2188   for (Int_t iter=0; iter<3; iter++){
2189     npoints=0;
2190     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2191       if (graph->GetY()[ipoint]==0) continue;
2192       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
2193       outx[npoints]  = graph->GetX()[ipoint];
2194       outy[npoints]  = graph->GetY()[ipoint];
2195       npoints++;
2196     }
2197     if (npoints<=1) break;
2198     medianY  =TMath::Median(npoints,outy);
2199     rmsY   =TMath::RMS(npoints,outy);
2200   }
2201   TGraph *graphOut=0;
2202   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2203   delete [] outx;
2204   delete [] outy;
2205   return graphOut;
2206 }
2207
2208
2209 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
2210   //
2211   // filter outlyer measurement
2212   // Only points around median +- cut filtered
2213   if (!graph) return  0;
2214   Int_t kMinPoints=2;
2215   Int_t npoints0 = graph->GetN();
2216   Int_t npoints=0;
2217   Float_t  rmsY=0;
2218   //
2219   //
2220   if (npoints0<kMinPoints) return 0;
2221
2222   Double_t *outx=new Double_t[npoints0];
2223   Double_t *outy=new Double_t[npoints0];
2224   for (Int_t iter=0; iter<3; iter++){
2225     npoints=0;
2226     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2227       if (graph->GetY()[ipoint]==0) continue;
2228       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
2229       outx[npoints]  = graph->GetX()[ipoint];
2230       outy[npoints]  = graph->GetY()[ipoint];
2231       npoints++;
2232     }
2233     if (npoints<=1) break;
2234     medianY  =TMath::Median(npoints,outy);
2235     rmsY   =TMath::RMS(npoints,outy);
2236   }
2237   TGraph *graphOut=0;
2238   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
2239   delete [] outx;
2240   delete [] outy;
2241   return graphOut;
2242 }
2243
2244
2245
2246 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * const graph, Float_t sigmaCut,Double_t &medianY){
2247   //
2248   // filter outlyer measurement
2249   // Only points with normalized errors median +- sigmaCut filtered
2250   //
2251   Int_t kMinPoints=10;
2252   Int_t npoints0 = graph->GetN();
2253   Int_t npoints=0;
2254   Float_t  medianErr=0, rmsErr=0;
2255   //
2256   //
2257   if (npoints0<kMinPoints) return 0;
2258
2259   Double_t *outx=new Double_t[npoints0];
2260   Double_t *outy=new Double_t[npoints0];
2261   Double_t *erry=new Double_t[npoints0];
2262   Double_t *nerry=new Double_t[npoints0];
2263   Double_t *errx=new Double_t[npoints0];
2264
2265   for (Int_t iter=0; iter<3; iter++){
2266     npoints=0;
2267     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
2268       nerry[npoints]  = graph->GetErrorY(ipoint);
2269       if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
2270       erry[npoints]  = graph->GetErrorY(ipoint);
2271       outx[npoints]  = graph->GetX()[ipoint];
2272       outy[npoints]  = graph->GetY()[ipoint];
2273       errx[npoints]  = graph->GetErrorY(ipoint);
2274       npoints++;
2275     }
2276     if (npoints==0) break;
2277     medianErr=TMath::Median(npoints,erry);
2278     medianY  =TMath::Median(npoints,outy);
2279     rmsErr   =TMath::RMS(npoints,erry);
2280   }
2281   TGraphErrors *graphOut=0;
2282   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
2283   delete []outx;
2284   delete []outy;
2285   delete []erry;
2286   delete []nerry;
2287   delete []errx;
2288   return graphOut;
2289 }
2290
2291
2292 void AliTPCcalibDButil::Sort(TGraph *graph){
2293   //
2294   // sort array - neccessay for approx
2295   //
2296   Int_t npoints = graph->GetN();
2297   Int_t *indexes=new Int_t[npoints];
2298   Double_t *outx=new Double_t[npoints];
2299   Double_t *outy=new Double_t[npoints];
2300   TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
2301   for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
2302   for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
2303   for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
2304   for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
2305
2306   delete [] indexes;
2307   delete [] outx;
2308   delete [] outy;
2309 }
2310 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
2311   //
2312   // smmoth graph - mean on the interval
2313   //
2314   Sort(graph);
2315   Int_t npoints = graph->GetN();
2316   Double_t *outy=new Double_t[npoints];
2317
2318   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2319     Double_t lx=graph->GetX()[ipoint];
2320     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2321     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2322     if (index0<0) index0=0;
2323     if (index1>=npoints-1) index1=npoints-1;
2324     if ((index1-index0)>1){
2325       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2326     }else{
2327       outy[ipoint]=graph->GetY()[ipoint];
2328     }
2329   }
2330  //  TLinearFitter  fitter(3,"pol2");
2331 //   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2332 //     Double_t lx=graph->GetX()[ipoint];
2333 //     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
2334 //     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
2335 //     if (index0<0) index0=0;
2336 //     if (index1>=npoints-1) index1=npoints-1;
2337 //     fitter.ClearPoints();
2338 //     for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
2339 //     if ((index1-index0)>1){
2340 //       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
2341 //     }else{
2342 //       outy[ipoint]=graph->GetY()[ipoint];
2343 //     }
2344 //   }
2345
2346
2347
2348   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
2349     graph->GetY()[ipoint] = outy[ipoint];
2350   }
2351   delete[] outy;
2352 }
2353
2354 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph * const graph, Double_t xref){
2355   //
2356   // Use constant interpolation outside of range
2357   //
2358   if (!graph) {
2359     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2360     return 0;
2361   }
2362
2363   if (graph->GetN()<1){
2364     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph \n");
2365     return 0;
2366   }
2367
2368
2369   if (xref<graph->GetX()[0]) return graph->GetY()[0];
2370   if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
2371
2372   //  printf("graph->Eval(graph->GetX()[0]) %f, graph->Eval(xref) %f \n",graph->Eval(graph->GetX()[0]), graph->Eval(xref));
2373
2374   if(graph->GetN()==1)
2375     return graph->Eval(graph->GetX()[0]);
2376
2377
2378   return graph->Eval(xref);
2379 }
2380
2381 Double_t AliTPCcalibDButil::EvalGraphConst(AliSplineFit *graph, Double_t xref){
2382   //
2383   // Use constant interpolation outside of range also for spline fits
2384   //
2385   if (!graph) {
2386     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
2387     return 0;
2388   }
2389   if (graph->GetKnots()<1){
2390     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
2391     return 0;
2392   }
2393   if (xref<graph->GetX()[0]) return graph->GetY0()[0];
2394   if (xref>graph->GetX()[graph->GetKnots()-1]) return graph->GetY0()[graph->GetKnots()-1];
2395   return graph->Eval( xref);
2396 }
2397
2398 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy,  Double_t sigmaCut){
2399   //
2400   // Filter DCS sensor information
2401   //   ymin     - minimal value
2402   //   ymax     - max value
2403   //   maxdy    - maximal deirivative
2404   //   sigmaCut - cut on values and derivative in terms of RMS distribution
2405   // Return value - accepted fraction
2406   //
2407   // Algorithm:
2408   //
2409   // 0. Calculate median and rms of values in specified range
2410   // 1. Filter out outliers - median+-sigmaCut*rms
2411   //    values replaced by median
2412   //
2413   AliSplineFit * fit    = sensor->GetFit();
2414   if (!fit) return 0.;
2415   Int_t          nknots = fit->GetKnots();
2416   if (nknots==0) {
2417     delete fit;
2418     sensor->SetFit(0);
2419     return 0;
2420   }
2421   //
2422   Double_t *yin0  = new Double_t[nknots];
2423   Double_t *yin1  = new Double_t[nknots];
2424   Int_t naccept=0;
2425
2426   for (Int_t iknot=0; iknot< nknots; iknot++){
2427     if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
2428       yin0[naccept]  = fit->GetY0()[iknot];
2429       yin1[naccept]  = fit->GetY1()[iknot];
2430       if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
2431       naccept++;
2432     }
2433   }
2434   if (naccept<1) {
2435     delete fit;
2436     sensor->SetFit(0);
2437     delete [] yin0;
2438     delete [] yin1;
2439     return 0.;
2440   }
2441
2442   Double_t medianY0=0, medianY1=0;
2443   Double_t rmsY0   =0, rmsY1=0;
2444   medianY0 = TMath::Median(naccept, yin0);
2445   medianY1 = TMath::Median(naccept, yin1);
2446   rmsY0    = TMath::RMS(naccept, yin0);
2447   rmsY1    = TMath::RMS(naccept, yin1);
2448   naccept=0;
2449   //
2450   // 1. Filter out outliers - median+-sigmaCut*rms
2451   //    values replaced by median
2452   //    if replaced the derivative set to 0
2453   //
2454   for (Int_t iknot=0; iknot< nknots; iknot++){
2455     Bool_t isOK=kTRUE;
2456     if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
2457     if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
2458     if (nknots<2) fit->GetY1()[iknot]=0;
2459     if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
2460     if (!isOK){
2461       fit->GetY0()[iknot]=medianY0;
2462       fit->GetY1()[iknot]=0;
2463     }else{
2464       naccept++;
2465     }
2466   }
2467   delete [] yin0;
2468   delete [] yin1;
2469   return Float_t(naccept)/Float_t(nknots);
2470 }
2471
2472 Float_t  AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
2473   //
2474   // Filter temperature array
2475   // tempArray    - array of temperatures         -
2476   // ymin         - minimal accepted temperature  - default 15
2477   // ymax         - maximal accepted temperature  - default 22
2478   // sigmaCut     - values filtered on interval median+-sigmaCut*rms - defaut 5
2479   // return value - fraction of filtered sensors
2480   const Double_t kMaxDy=0.1;
2481   Int_t nsensors=tempArray->NumSensors();
2482   if (nsensors==0) return 0.;
2483   Int_t naccept=0;
2484   for (Int_t isensor=0; isensor<nsensors; isensor++){
2485     AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
2486     if (!sensor) continue;
2487     //printf("%d\n",isensor);
2488     FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
2489     if (sensor->GetFit()==0){
2490       //delete sensor;
2491       tempArray->RemoveSensorNum(isensor);
2492     }else{
2493       naccept++;
2494     }
2495   }
2496   return Float_t(naccept)/Float_t(nsensors);
2497 }
2498
2499
2500 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector * const pcstream){
2501   //
2502   // Filter CE data
2503   // Input parameters:
2504   //    deltaT   - smoothing window (in seconds)
2505   //    cutAbs   - max distance of the time info to the median (in time bins)
2506   //    cutSigma - max distance (in the RMS)
2507   //    pcstream - optional debug streamer to store original and filtered info
2508   // Hardwired parameters:
2509   //    kMinPoints =10;       // minimal number of points to define the CE
2510   //    kMinSectors=12;       // minimal number of sectors to define sideCE
2511   // Algorithm:
2512   // 0. Filter almost emty graphs (kMinPoints=10)
2513   // 1. calculate median and RMS per side
2514   // 2. Filter graphs - in respect with side medians
2515   //                  - cutAbs and cutDelta used
2516   // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
2517   // 4. Calculate mean for A side and C side
2518   //
2519   const Int_t kMinPoints =10;       // minimal number of points to define the CE
2520   const Int_t kMinSectors=12;       // minimal number of sectors to define sideCE
2521   const Int_t kMinTime   =400;     // minimal arrival time of CE
2522   TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
2523   Double_t medianY=0;
2524   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
2525   if (!cearray) return;
2526   Double_t tmin=-1;
2527   Double_t tmax=-1;
2528   //
2529   //
2530   AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
2531   AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernAtmosPressure");
2532   if ( tempMapCE && cavernPressureCE){
2533     //
2534     //     Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
2535     //     FilterSensor(cavernPressureCE,960,1050,10, 5.);
2536     //     if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
2537     Bool_t isOK=kTRUE;
2538     if (isOK)  {
2539       // recalculate P/T correction map for time of the CE
2540       AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
2541       driftCalib->SetName("driftPTCE");
2542       driftCalib->SetTitle("driftPTCE");
2543       cearray->AddLast(driftCalib);
2544     }
2545   }
2546   //
2547   // 0. Filter almost emty graphs
2548   //
2549
2550   for (Int_t i=0; i<72;i++){
2551     TGraph *graph= (TGraph*)arrT->At(i);
2552     if (!graph) continue;
2553     graph->Sort();
2554     if (graph->GetN()<kMinPoints){
2555       arrT->AddAt(0,i);
2556       delete graph;  // delete empty graph
2557       continue;
2558     }
2559     if (tmin<0) tmin = graph->GetX()[0];
2560     if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
2561     //
2562     if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
2563     if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
2564   }
2565   //
2566   // 1. calculate median and RMS per side
2567   //
2568   TArrayF arrA(100000), arrC(100000);
2569   Int_t nA=0, nC=0;
2570   Double_t medianA=0, medianC=0;
2571   Double_t rmsA=0, rmsC=0;
2572   for (Int_t isec=0; isec<72;isec++){
2573     TGraph *graph= (TGraph*)arrT->At(isec);
2574     if (!graph) continue;
2575     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
2576       if (graph->GetY()[ipoint]<kMinTime) continue;
2577       if (nA>=arrA.fN) arrA.Set(nA*2);
2578       if (nC>=arrC.fN) arrC.Set(nC*2);
2579       if (isec%36<18)  arrA[nA++]= graph->GetY()[ipoint];
2580       if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
2581     }
2582   }
2583   if (nA>0){
2584     medianA=TMath::Median(nA,arrA.fArray);
2585     rmsA   =TMath::RMS(nA,arrA.fArray);
2586   }
2587   if (nC>0){
2588     medianC=TMath::Median(nC,arrC.fArray);
2589     rmsC   =TMath::RMS(nC,arrC.fArray);
2590   }
2591   //
2592   // 2. Filter graphs - in respect with side medians
2593   //
2594   TArrayD vecX(100000), vecY(100000);
2595   for (Int_t isec=0; isec<72;isec++){
2596     TGraph *graph= (TGraph*)arrT->At(isec);
2597     if (!graph) continue;
2598     Double_t median = (isec%36<18) ? medianA: medianC;
2599     Double_t rms    = (isec%36<18) ? rmsA:    rmsC;
2600     Int_t naccept=0;
2601     //    for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){ //not neccessary to remove first points
2602     for (Int_t ipoint=0; ipoint<graph->GetN();ipoint++){
2603       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
2604       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
2605       vecX[naccept]= graph->GetX()[ipoint];
2606       vecY[naccept]= graph->GetY()[ipoint];
2607       naccept++;
2608     }
2609     if (naccept<kMinPoints){
2610       arrT->AddAt(0,isec);
2611       delete graph;  // delete empty graph
2612       continue;
2613     }
2614     TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
2615     delete graph;
2616     arrT->AddAt(graph2,isec);
2617   }
2618   //
2619   // 3. Cut in respect wit the graph median
2620   //
2621   for (Int_t i=0; i<72;i++){
2622     TGraph *graph= (TGraph*)arrT->At(i);
2623     if (!graph) continue;
2624     //
2625     // filter in range
2626     //
2627     TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2628     if (!graphTS0) continue;
2629     if (graphTS0->GetN()<kMinPoints) {
2630       delete graphTS0;
2631       delete graph;
2632       arrT->AddAt(0,i);
2633       continue;
2634     }
2635     TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2636     graphTS->Sort();
2637     AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2638     if (pcstream){
2639       Int_t run = AliTPCcalibDB::Instance()->GetRun();
2640       (*pcstream)<<"filterCE"<<
2641         "run="<<run<<
2642         "isec="<<i<<
2643         "mY="<<medianY<<
2644         "graph.="<<graph<<
2645         "graphTS0.="<<graphTS0<<
2646         "graphTS.="<<graphTS<<
2647         "\n";
2648     }
2649     delete graphTS0;
2650     if (!graphTS) continue;
2651     arrT->AddAt(graphTS,i);
2652     delete graph;
2653   }
2654   //
2655   // Recalculate the mean time A side C side
2656   //
2657   TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2658   Int_t meanPoints=(nA+nC)/72;  // mean number of points
2659   for (Int_t itime=0; itime<200; itime++){
2660     nA=0, nC=0;
2661     Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2662     for (Int_t i=0; i<72;i++){
2663       TGraph *graph= (TGraph*)arrT->At(i);
2664       if (!graph) continue;
2665       if (graph->GetN()<(meanPoints/4)) continue;
2666       if ( (i%36)<18 )  arrA[nA++]=graph->Eval(time);
2667       if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2668     }
2669     xA[itime]=time;
2670     xC[itime]=time;
2671     yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2672     yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2673     eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2674     eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2675   }
2676   //
2677   Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2678   Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2679   if (pcstream){
2680     Int_t run = AliTPCcalibDB::Instance()->GetRun();
2681     (*pcstream)<<"filterAC"<<
2682       "run="<<run<<
2683       "nA="<<nA<<
2684       "nC="<<nC<<
2685       "rmsTA="<<rmsTA<<
2686       "rmsTC="<<rmsTC<<
2687       "\n";
2688   }
2689   //
2690   TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2691   TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2692   TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2693   if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2694   TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2695   if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2696   delete grA;
2697   delete grC;
2698   if (nA<kMinSectors) arrT->AddAt(0,72);
2699   else arrT->AddAt(graphTSA,72);
2700   if (nC<kMinSectors) arrT->AddAt(0,73);
2701   else arrT->AddAt(graphTSC,73);
2702 }
2703
2704
2705 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector * const pcstream){
2706   //
2707   // Filter Drift velocity measurement using the tracks
2708   // 0.  remove outlyers - error based
2709   //     cutSigma
2710   //
2711   //
2712   const Int_t kMinPoints=1;  // minimal number of points to define value
2713   TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2714   Double_t medianY=0;
2715   if (!arrT) return;
2716   for (Int_t i=0; i<arrT->GetEntries();i++){
2717     TGraphErrors *graph= dynamic_cast<TGraphErrors*>(arrT->At(i));
2718     if (!graph) continue;
2719     if (graph->GetN()<kMinPoints){
2720       delete graph;
2721       arrT->AddAt(0,i);
2722       continue;
2723     }
2724     TGraphErrors *graph2 = NULL;
2725     if(graph->GetN()<10) {
2726       graph2 = new TGraphErrors(graph->GetN(),graph->GetX(),graph->GetY(),graph->GetEX(),graph->GetEY());
2727       if (!graph2) {
2728         delete graph; arrT->AddAt(0,i); continue;
2729       }
2730     }
2731     else {
2732       graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2733       if (!graph2) {
2734         delete graph; arrT->AddAt(0,i); continue;
2735       }
2736     }
2737     if (graph2->GetN()<1) {
2738       delete graph; arrT->AddAt(0,i); continue;
2739     }
2740     graph2->SetName(graph->GetName());
2741     graph2->SetTitle(graph->GetTitle());
2742     arrT->AddAt(graph2,i);
2743     if (pcstream){
2744       (*pcstream)<<"filterTracks"<<
2745         "run="<<run<<
2746         "isec="<<i<<
2747         "mY="<<medianY<<
2748         "graph.="<<graph<<
2749         "graph2.="<<graph2<<
2750         "\n";
2751     }
2752     delete graph;
2753   }
2754 }
2755
2756
2757
2758
2759
2760 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2761   //
2762   //
2763   // get laser time offset
2764   // median around timeStamp+-deltaT
2765   // QA - chi2 needed for later usage - to be added
2766   //    - currently cut on error
2767   //
2768   Int_t kMinPoints=1;
2769   Double_t kMinDelay=0.01;
2770   Double_t kMinDelayErr=0.0001;
2771
2772   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2773   if (!array) return 0;
2774   TGraphErrors *tlaser=0;
2775   if (array){
2776     if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2777     if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2778   }
2779   if (!tlaser) return 0;
2780   Int_t npoints0= tlaser->GetN();
2781   if (npoints0==0) return 0;
2782   Double_t *xlaser = new Double_t[npoints0];
2783   Double_t *ylaser = new Double_t[npoints0];
2784   Int_t npoints=0;
2785   for (Int_t i=0;i<npoints0;i++){
2786     //printf("%d\n",i);
2787     if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2788     if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2789     xlaser[npoints]=tlaser->GetX()[npoints];
2790     ylaser[npoints]=tlaser->GetY()[npoints];
2791     npoints++;
2792   }
2793   //
2794   //
2795   Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2796   Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2797   //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2798   if (index0<0) index0=0;
2799   if (index1>=npoints-1) index1=npoints-1;
2800   if (index1-index0<kMinPoints) return 0;
2801   //
2802   //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2803     Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2804   delete [] ylaser;
2805   delete [] xlaser;
2806   return mean;
2807 }
2808
2809
2810
2811
2812 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector * const pcstream){
2813   //
2814   // Filter Goofie data
2815   // goofieArray - points will be filtered
2816   // deltaT      - smmothing time window
2817   // cutSigma    - outler sigma cut in rms
2818   // minVn, maxVd- range absolute cut for variable vd/pt
2819   //             - to be tuned
2820   //
2821   // Ignore goofie if not enough points
2822   //
2823   const Int_t kMinPoints = 3;
2824   //
2825
2826   TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2827   TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2828   TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2829   TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2830   if (!graphvd) return;
2831   if (graphvd->GetN()<kMinPoints){
2832     delete graphvd;
2833     goofieArray->GetSensorNum(2)->SetGraph(0);
2834     return;
2835   }
2836   //
2837   // 1. Caluclate medians of critical variables
2838   //    drift velcocity
2839   //    P/T
2840   //    area near peak
2841   //    area far  peak
2842   //
2843   Double_t medianpt=0;
2844   Double_t medianvd=0, sigmavd=0;
2845   Double_t medianan=0;
2846   Double_t medianaf=0;
2847   Int_t    entries=graphvd->GetN();
2848   Double_t yvdn[10000];
2849   Int_t nvd=0;
2850   //
2851   for (Int_t ipoint=0; ipoint<entries; ipoint++){
2852     if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2853     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2854     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2855     yvdn[nvd++]=graphvd->GetY()[ipoint];
2856   }
2857   if (nvd<kMinPoints){
2858     delete graphvd;
2859     goofieArray->GetSensorNum(2)->SetGraph(0);
2860     return;
2861   }
2862   //
2863   Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2864   if (nuni>=kMinPoints){
2865     AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2866   }else{
2867     medianvd = TMath::Median(nvd, yvdn);
2868   }
2869
2870   TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2871   TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2872   TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2873   TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2874   TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2875   TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2876   delete graphpt0;
2877   delete graphpt1;
2878   delete graphan0;
2879   delete graphan1;
2880   delete graphaf0;
2881   delete graphaf1;
2882   //
2883   // 2. Make outlyer graph
2884   //
2885   Int_t nOK=0;
2886   TGraph graphOut(*graphvd);
2887   for (Int_t i=0; i<entries;i++){
2888     //
2889     Bool_t isOut=kFALSE;
2890     if (graphpt->GetY()[i]<=0.0000001) {  graphOut.GetY()[i]=1; continue;}
2891     if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) {  graphOut.GetY()[i]=1; continue;}
2892
2893     if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2894       isOut|=kTRUE;
2895     if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2896     if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2897     if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2898     graphOut.GetY()[i]= (isOut)?1:0;
2899     if (!isOut) nOK++;
2900   }
2901   if (nOK<kMinPoints) {
2902     delete graphvd;
2903     goofieArray->GetSensorNum(2)->SetGraph(0);
2904     return;
2905   }
2906   //
2907   // 3. Filter out outlyers - and smooth
2908   //
2909   TVectorF vmedianArray(goofieArray->NumSensors());
2910   TVectorF vrmsArray(goofieArray->NumSensors());
2911   Double_t xnew[10000];
2912   Double_t ynew[10000];
2913   TObjArray junk;
2914   junk.SetOwner(kTRUE);
2915   Bool_t isOK=kTRUE;
2916   //
2917   //
2918   for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2919     isOK=kTRUE;
2920     AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2921     TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2922     //
2923     if (!sensor) continue;
2924     graphOld = sensor->GetGraph();
2925     if (graphOld) {
2926       sensor->SetGraph(0);
2927       Int_t nused=0;
2928       for (Int_t i=0;i<entries;i++){
2929         if (graphOut.GetY()[i]>0.5) continue;
2930         xnew[nused]=graphOld->GetX()[i];
2931         ynew[nused]=graphOld->GetY()[i];
2932         nused++;
2933       }
2934       graphNew = new TGraph(nused,xnew,ynew);
2935       junk.AddLast(graphNew);
2936       junk.AddLast(graphOld);
2937       Double_t median=0;
2938       graphNew0  = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2939       if (graphNew0!=0){
2940         junk.AddLast(graphNew0);
2941         graphNew1  = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2942         if (graphNew1!=0){
2943           junk.AddLast(graphNew1);
2944           graphNew2  = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2945           if (graphNew2!=0) {
2946             vrmsArray[isensor]   =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2947             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2948             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2949             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2950             printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2951             vmedianArray[isensor]=median;
2952             //
2953           }
2954         }
2955       }
2956     }
2957     if (!graphOld) {  isOK=kFALSE; graphOld =&graphOut;}
2958     if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2959     if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2960     if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2961     (*pcstream)<<"goofieA"<<
2962       Form("isOK_%d.=",isensor)<<isOK<<
2963       Form("s_%d.=",isensor)<<sensor<<
2964       Form("gr_%d.=",isensor)<<graphOld<<
2965       Form("gr0_%d.=",isensor)<<graphNew0<<
2966       Form("gr1_%d.=",isensor)<<graphNew1<<
2967       Form("gr2_%d.=",isensor)<<graphNew2;
2968     if (isOK) sensor->SetGraph(graphNew2);
2969   }
2970   (*pcstream)<<"goofieA"<<
2971     "vmed.="<<&vmedianArray<<
2972     "vrms.="<<&vrmsArray<<
2973     "\n";
2974   junk.Delete();   // delete temoprary graphs
2975
2976 }
2977
2978
2979
2980
2981
2982 TMatrixD* AliTPCcalibDButil::MakeStatRelKalman(TObjArray * const array, Float_t minFraction, Int_t minStat, Float_t maxvd){
2983   //
2984   // Make a statistic matrix
2985   // Input parameters:
2986   //   array        - TObjArray of AliRelKalmanAlign
2987   //   minFraction  - minimal ration of accepted tracks
2988   //   minStat      - minimal statistic (number of accepted tracks)
2989   //   maxvd        - maximal deviation for the 1
2990   // Output matrix:
2991   //    columns    - Mean, Median, RMS
2992   //    row        - parameter type (rotation[3], translation[3], drift[3])
2993   if (!array) return 0;
2994   if (array->GetEntries()<=0) return 0;
2995   //  Int_t entries = array->GetEntries();
2996   Int_t entriesFast = array->GetEntriesFast();
2997   TVectorD state(9);
2998   TVectorD *valArray[9];
2999   for (Int_t i=0; i<9; i++){
3000     valArray[i] = new TVectorD(entriesFast);
3001   }
3002   Int_t naccept=0;
3003   for (Int_t ikalman=0; ikalman<entriesFast; ikalman++){
3004     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(ikalman);
3005     if (!kalman) continue;
3006     if (TMath::Abs(kalman->GetTPCvdCorr()-1)>maxvd) continue;
3007     if (kalman->GetNUpdates()<minStat) continue;
3008     if (kalman->GetNUpdates()/kalman->GetNTracks()<minFraction) continue;
3009     kalman->GetState(state);
3010     for (Int_t ipar=0; ipar<9; ipar++)
3011       (*valArray[ipar])[naccept]=state[ipar];
3012     naccept++;
3013   }
3014   //if (naccept<2) return 0;
3015   if (naccept<1) return 0;
3016   TMatrixD *pstat=new TMatrixD(9,3);
3017   TMatrixD &stat=*pstat;
3018   for (Int_t ipar=0; ipar<9; ipar++){
3019     stat(ipar,0)=TMath::Mean(naccept, valArray[ipar]->GetMatrixArray());
3020     stat(ipar,1)=TMath::Median(naccept, valArray[ipar]->GetMatrixArray());
3021     stat(ipar,2)=TMath::RMS(naccept, valArray[ipar]->GetMatrixArray());
3022   }
3023   return pstat;
3024 }
3025
3026
3027 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const array, const TMatrixD & stat, Bool_t direction, Float_t sigmaCut){
3028   //
3029   // Smooth the array of AliRelKalmanAlign - detector alignment and drift calibration)
3030   // Input:
3031   //   array     - input array
3032   //   stat      - mean parameters statistic
3033   //   direction -
3034   //   sigmaCut  - maximal allowed deviation from mean in terms of RMS
3035   if (!array) return 0;
3036   if (array->GetEntries()<=0) return 0;
3037   if (!(&stat)) return 0;
3038   // error increase in 1 hour
3039   const Double_t kerrsTime[9]={
3040     0.00001,  0.00001, 0.00001,
3041     0.001,    0.001,   0.001,
3042     0.002,  0.01,   0.001};
3043   //
3044   //
3045   Int_t entries = array->GetEntriesFast();
3046   TObjArray *sArray= new TObjArray(entries);
3047   AliRelAlignerKalman * sKalman =0;
3048   TVectorD state(9);
3049   for (Int_t i=0; i<entries; i++){
3050     Int_t index=(direction)? entries-i-1:i;
3051     AliRelAlignerKalman * kalman = (AliRelAlignerKalman *) array->UncheckedAt(index);
3052     if (!kalman) continue;
3053     Bool_t isOK=kTRUE;
3054     kalman->GetState(state);
3055     for (Int_t ipar=0; ipar<9; ipar++){
3056       if (TMath::Abs(state[ipar]-stat(ipar,1))>sigmaCut*stat(ipar,2)) isOK=kFALSE;
3057     }
3058     if (!sKalman &&isOK) {
3059       sKalman=new AliRelAlignerKalman(*kalman);
3060       sKalman->SetRejectOutliers(kFALSE);
3061       sKalman->SetRunNumber(kalman->GetRunNumber());
3062       sKalman->SetTimeStamp(kalman->GetTimeStamp());
3063     }
3064     if (!sKalman) continue;
3065     Double_t deltaT=TMath::Abs(Int_t(kalman->GetTimeStamp())-Int_t(sKalman->GetTimeStamp()))/3600.;
3066     for (Int_t ipar=0; ipar<9; ipar++){
3067 //       (*(sKalman->GetStateCov()))(6,6)+=deltaT*errvd*errvd;
3068 //       (*(sKalman->GetStateCov()))(7,7)+=deltaT*errt0*errt0;
3069 //       (*(sKalman->GetStateCov()))(8,8)+=deltaT*errvy*errvy;
3070       (*(sKalman->GetStateCov()))(ipar,ipar)+=deltaT*kerrsTime[ipar]*kerrsTime[ipar];
3071     }
3072     sKalman->SetRunNumber(kalman->GetRunNumber());
3073     if (!isOK) sKalman->SetRunNumber(0);
3074     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3075     if (!isOK) continue;
3076     sKalman->SetRejectOutliers(kFALSE);
3077     sKalman->SetRunNumber(kalman->GetRunNumber());
3078     sKalman->SetTimeStamp(kalman->GetTimeStamp());
3079     sKalman->Merge(kalman);
3080     sArray->AddAt(new AliRelAlignerKalman(*sKalman),index);
3081     //sKalman->Print();
3082   }
3083   return sArray;
3084 }
3085
3086 TObjArray *AliTPCcalibDButil::SmoothRelKalman(TObjArray * const arrayP, TObjArray * const arrayM){
3087   //
3088   // Merge 2 RelKalman arrays
3089   // Input:
3090   //   arrayP    - rel kalman in direction plus
3091   //   arrayM    - rel kalman in direction minus
3092   if (!arrayP) return 0;
3093   if (arrayP->GetEntries()<=0) return 0;
3094   if (!arrayM) return 0;
3095   if (arrayM->GetEntries()<=0) return 0;
3096
3097   Int_t entries = arrayP->GetEntriesFast();
3098   TObjArray *array = new TObjArray(arrayP->GetEntriesFast());
3099
3100   for (Int_t i=0; i<entries; i++){
3101      AliRelAlignerKalman * kalmanP = (AliRelAlignerKalman *) arrayP->UncheckedAt(i);
3102      AliRelAlignerKalman * kalmanM = (AliRelAlignerKalman *) arrayM->UncheckedAt(i);
3103      if (!kalmanP) continue;
3104      if (!kalmanM) continue;
3105
3106      AliRelAlignerKalman *kalman = NULL;
3107      if(kalmanP->GetRunNumber() != 0 && kalmanM->GetRunNumber() != 0) {
3108        kalman = new AliRelAlignerKalman(*kalmanP);
3109        kalman->Merge(kalmanM);
3110      }
3111      else if (kalmanP->GetRunNumber() == 0) {
3112        kalman = new AliRelAlignerKalman(*kalmanM);
3113      }
3114      else if (kalmanM->GetRunNumber() == 0) {
3115        kalman = new AliRelAlignerKalman(*kalmanP);
3116      }
3117      else
3118        continue;
3119
3120      array->AddAt(kalman,i);
3121   }
3122
3123   return array;
3124 }