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