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