TPC/AliTPCcalibDButil.cxx

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