TPC/AliTPCcalibDButil.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // Class providing the calculation of derived quantities (mean,rms,fits,...) //
  20 //       of calibration entries                                              //
  21 /*
  22
  23
  24 */
  25 ////////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TMath.h>
  28 #include <TVectorT.h>
  29 #include <TObjArray.h>
  30 #include <TGraph.h>
  31 #include <TFile.h>
  32 #include <TDirectory.h>
  33 #include <AliCDBStorage.h>
  34 #include <AliDCSSensorArray.h>
  35 #include <AliTPCSensorTempArray.h>
  36 #include <AliDCSSensor.h>
  37 #include "AliTPCcalibDB.h"
  38 #include "AliTPCCalPad.h"
  39 #include "AliTPCCalROC.h"
  40 #include "AliTPCROC.h"
  41 #include "AliTPCmapper.h"
  42 #include "AliTPCParam.h"
  43 #include "AliTPCCalibRaw.h"
  44 #include "TGraphErrors.h"
  45 #include "AliLog.h"
  46 #include "AliTPCcalibDButil.h"
  47 #include "AliTPCPreprocessorOnline.h"
  48 #include "AliTPCCalibVdrift.h"
  49 #include "AliMathBase.h"
  50
  51 ClassImp(AliTPCcalibDButil)
  52 AliTPCcalibDButil::AliTPCcalibDButil() :
  53   TObject(),
  54   fCalibDB(0),
  55   fPadNoise(0x0),
  56   fPedestals(0x0),
  57   fPulserTmean(0x0),
  58   fPulserTrms(0x0),
  59   fPulserQmean(0x0),
  60   fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
  61   fCETmean(0x0),
  62   fCETrms(0x0),
  63   fCEQmean(0x0),
  64   fALTROMasked(0x0),
  65   fCalibRaw(0x0),
  66   fRefPadNoise(0x0),
  67   fRefPedestals(0x0),
  68   fRefPulserTmean(0x0),
  69   fRefPulserTrms(0x0),
  70   fRefPulserQmean(0x0),
  71   fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
  72   fRefCETmean(0x0),
  73   fRefCETrms(0x0),
  74   fRefCEQmean(0x0),
  75   fRefALTROMasked(0x0),
  76   fRefCalibRaw(0x0),
  77   fGoofieArray(0x0),
  78   fMapper(new AliTPCmapper(0x0)),
  79   fNpulserOutliers(-1),
  80   fIrocTimeOffset(0),
  81   fCETmaxLimitAbs(1.5),
  82   fPulTmaxLimitAbs(1.5),
  83   fPulQmaxLimitAbs(5),
  84   fPulQminLimit(11),
  85   fRuns(0),                         // run list with OCDB info
  86   fRunsStart(0),                    // start time for given run
  87   fRunsStop(0)                     // stop time for given run
  88 {
  89   //
  90   // Default ctor
  91   //
  92 }
  93 //_____________________________________________________________________________________
  94 AliTPCcalibDButil::~AliTPCcalibDButil()
  95 {
  96   //
  97   // dtor
  98   //
  99   delete fPulserOutlier;
 100   delete fRefPulserOutlier;
 101   delete fMapper;
 102   if (fRefPadNoise) delete fRefPadNoise;
 103   if (fRefPedestals) delete fRefPedestals;
 104   if (fRefPulserTmean) delete fRefPulserTmean;
 105   if (fRefPulserTrms) delete fRefPulserTrms;
 106   if (fRefPulserQmean) delete fRefPulserQmean;
 107   if (fRefCETmean) delete fRefCETmean;
 108   if (fRefCETrms) delete fRefCETrms;
 109   if (fRefCEQmean) delete fRefCEQmean;
 110   if (fRefALTROMasked) delete fRefALTROMasked;
 111   if (fRefCalibRaw) delete fRefCalibRaw;
 112
 113 }
 114 //_____________________________________________________________________________________
 115 void AliTPCcalibDButil::UpdateFromCalibDB()
 116 {
 117   //
 118   // Update pointers from calibDB
 119   //
 120   if (!fCalibDB) fCalibDB=AliTPCcalibDB::Instance();
 121   fPadNoise=fCalibDB->GetPadNoise();
 122   fPedestals=fCalibDB->GetPedestals();
 123   fPulserTmean=fCalibDB->GetPulserTmean();
 124   fPulserTrms=fCalibDB->GetPulserTrms();
 125   fPulserQmean=fCalibDB->GetPulserQmean();
 126   fCETmean=fCalibDB->GetCETmean();
 127   fCETrms=fCalibDB->GetCETrms();
 128   fCEQmean=fCalibDB->GetCEQmean();
 129   fALTROMasked=fCalibDB->GetALTROMasked();
 130   fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
 131   fCalibRaw=fCalibDB->GetCalibRaw();
 132   UpdatePulserOutlierMap();
 133 }
 134 //_____________________________________________________________________________________
 135 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
 136                                       Int_t &noutliersCE, Double_t & chi2A, Double_t &chi2C, AliTPCCalPad *outCE)
 137 {
 138   //
 139   // Process the CE data for this run
 140   // the return TVectorD arrays contian the results of the fit
 141   // noutliersCE contains the number of pads marked as outliers,
 142   //   not including masked and edge pads
 143   //
 144
 145   //retrieve CE and ALTRO data
 146   if (!fCETmean){
 147     TString fitString(fitFormula);
 148     fitString.ReplaceAll("++","#");
 149     Int_t ndim=fitString.CountChar('#')+2;
 150     fitResultsA.ResizeTo(ndim);
 151     fitResultsC.ResizeTo(ndim);
 152     fitResultsA.Zero();
 153     fitResultsC.Zero();
 154     noutliersCE=-1;
 155     return;
 156   }
 157   noutliersCE=0;
 158   //create outlier map
 159   AliTPCCalPad *out=0;
 160   if (outCE) out=outCE;
 161   else out=new AliTPCCalPad("outCE","outCE");
 162   AliTPCCalROC *rocMasked=0x0;
 163   //loop over all channels
 164   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 165     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 166     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
 167     AliTPCCalROC *rocOut=out->GetCalROC(iroc);
 168     if (!rocData) {
 169       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
 170       rocOut->Add(1.);
 171       continue;
 172     }
 173     //add time offset to IROCs
 174     if (iroc<AliTPCROC::Instance()->GetNInnerSector())
 175       rocData->Add(fIrocTimeOffset);
 176     //select outliers
 177     UInt_t nrows=rocData->GetNrows();
 178     for (UInt_t irow=0;irow<nrows;++irow){
 179       UInt_t npads=rocData->GetNPads(irow);
 180       for (UInt_t ipad=0;ipad<npads;++ipad){
 181         rocOut->SetValue(irow,ipad,0);
 182         //exclude masked pads
 183         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
 184           rocOut->SetValue(irow,ipad,1);
 185           continue;
 186         }
 187         //exclude first two rows in IROC and last two rows in OROC
 188         if (iroc<36){
 189           if (irow<2) rocOut->SetValue(irow,ipad,1);
 190         } else {
 191           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
 192         }
 193         //exclude edge pads
 194         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
 195         Float_t valTmean=rocData->GetValue(irow,ipad);
 196         //exclude values that are exactly 0
 197         if (valTmean==0) {
 198           rocOut->SetValue(irow,ipad,1);
 199           ++noutliersCE;
 200         }
 201         // exclude channels with too large variations
 202         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
 203           rocOut->SetValue(irow,ipad,1);
 204           ++noutliersCE;
 205         }
 206       }
 207     }
 208   }
 209   //perform fit
 210   TMatrixD dummy;
 211   Float_t chi2Af,chi2Cf;
 212   fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
 213   chi2A=chi2Af;
 214   chi2C=chi2Cf;
 215   if (!outCE) delete out;
 216 }
 217 //_____________________________________________________________________________________
 218 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
 219                      TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
 220                      Float_t &driftTimeA, Float_t &driftTimeC )
 221 {
 222   //
 223   // Calculate statistical information from the CE graphs for drift time and charge
 224   //
 225
 226   //reset arrays
 227   vecTEntries.ResizeTo(72);
 228   vecTMean.ResizeTo(72);
 229   vecTRMS.ResizeTo(72);
 230   vecTMedian.ResizeTo(72);
 231   vecQEntries.ResizeTo(72);
 232   vecQMean.ResizeTo(72);
 233   vecQRMS.ResizeTo(72);
 234   vecQMedian.ResizeTo(72);
 235   vecTEntries.Zero();
 236   vecTMean.Zero();
 237   vecTRMS.Zero();
 238   vecTMedian.Zero();
 239   vecQEntries.Zero();
 240   vecQMean.Zero();
 241   vecQRMS.Zero();
 242   vecQMedian.Zero();
 243   driftTimeA=0;
 244   driftTimeC=0;
 245   TObjArray *arrT=fCalibDB->GetCErocTtime();
 246   TObjArray *arrQ=fCalibDB->GetCErocQtime();
 247   if (arrT){
 248     for (Int_t isec=0;isec<74;++isec){
 249       TGraph *gr=(TGraph*)arrT->At(isec);
 250       if (!gr) continue;
 251       TVectorD values;
 252       Int_t npoints = gr->GetN();
 253       values.ResizeTo(npoints);
 254       Int_t nused =0;
 255       //skip first points, theres always some problems with finding the CE position
 256       for (Int_t ipoint=4; ipoint<npoints; ipoint++){
 257         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
 258           values[nused]=gr->GetY()[ipoint];
 259           nused++;
 260         }
 261       }
 262       //
 263       if (isec<72) vecTEntries[isec]= nused;
 264       if (nused>1){
 265         if (isec<72){
 266           vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 267           vecTMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 268           vecTRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 269         } else if (isec==72){
 270           driftTimeA=TMath::Median(nused,values.GetMatrixArray());
 271         } else if (isec==73){
 272           driftTimeC=TMath::Median(nused,values.GetMatrixArray());
 273         }
 274       }
 275     }
 276   }
 277   if (arrQ){
 278     for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
 279       TGraph *gr=(TGraph*)arrQ->At(isec);
 280       if (!gr) continue;
 281       TVectorD values;
 282       Int_t npoints = gr->GetN();
 283       values.ResizeTo(npoints);
 284       Int_t nused =0;
 285       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 286         if (gr->GetY()[ipoint]>10 && gr->GetY()[ipoint]<500 ){
 287           values[nused]=gr->GetY()[ipoint];
 288           nused++;
 289         }
 290       }
 291       //
 292       vecQEntries[isec]= nused;
 293       if (nused>1){
 294         vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 295         vecQMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 296         vecQRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 297       }
 298     }
 299   }
 300 }
 301
 302 //_____________________________________________________________________________________
 303 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
 304                       TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
 305                       Int_t &nonMaskedZero)
 306 {
 307   //
 308   // process noise data
 309   // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
 310   //    OROCs small pads [2] and OROCs large pads [3]
 311   // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
 312   // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
 313   //
 314
 315   //set proper size and reset
 316   const UInt_t infoSize=4;
 317   vNoiseMean.ResizeTo(infoSize);
 318   vNoiseMeanSenRegions.ResizeTo(infoSize);
 319   vNoiseRMS.ResizeTo(infoSize);
 320   vNoiseRMSSenRegions.ResizeTo(infoSize);
 321   vNoiseMean.Zero();
 322   vNoiseMeanSenRegions.Zero();
 323   vNoiseRMS.Zero();
 324   vNoiseRMSSenRegions.Zero();
 325   nonMaskedZero=0;
 326   //counters
 327   TVectorD c(infoSize);
 328   TVectorD cs(infoSize);
 329   //tpc parameters
 330   AliTPCParam par;
 331   par.Update();
 332   //retrieve noise and ALTRO data
 333   if (!fPadNoise) return;
 334   AliTPCCalROC *rocMasked=0x0;
 335   //create IROC, OROC1, OROC2 and sensitive region masks
 336   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 337     AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
 338     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 339     UInt_t nrows=noiseROC->GetNrows();
 340     for (UInt_t irow=0;irow<nrows;++irow){
 341       UInt_t npads=noiseROC->GetNPads(irow);
 342       for (UInt_t ipad=0;ipad<npads;++ipad){
 343         //don't use masked channels;
 344         if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
 345         Float_t noiseVal=noiseROC->GetValue(irow,ipad);
 346         //check if noise==0
 347         if (noiseVal==0) {
 348           ++nonMaskedZero;
 349           continue;
 350         }
 351         //check for nan
 352         if ( !(noiseVal<10000000) ){
 353           printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
 354           continue;
 355         }
 356         Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
 357         Int_t masksen=1; // sensitive pards are not masked (0)
 358         if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
 359         if (isec<AliTPCROC::Instance()->GetNInnerSector()){
 360           //IROCs
 361           if (irow>19&&irow<46){
 362             if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
 363           }
 364           Int_t type=1;
 365           vNoiseMean[type]+=noiseVal;
 366           vNoiseRMS[type]+=noiseVal*noiseVal;
 367           ++c[type];
 368           if (!masksen){
 369             vNoiseMeanSenRegions[type]+=noiseVal;
 370             vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 371             ++cs[type];
 372           }
 373         } else {
 374           //OROCs
 375           //define sensive regions
 376           if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
 377           if ( irow>75 ){
 378             Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
 379             if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
 380           }
 381           if ((Int_t)irow<par.GetNRowUp1()){
 382             //OROC1
 383             Int_t type=2;
 384             vNoiseMean[type]+=noiseVal;
 385             vNoiseRMS[type]+=noiseVal*noiseVal;
 386             ++c[type];
 387             if (!masksen){
 388               vNoiseMeanSenRegions[type]+=noiseVal;
 389               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 390               ++cs[type];
 391             }
 392           }else{
 393             //OROC2
 394             Int_t type=3;
 395             vNoiseMean[type]+=noiseVal;
 396             vNoiseRMS[type]+=noiseVal*noiseVal;
 397             ++c[type];
 398             if (!masksen){
 399               vNoiseMeanSenRegions[type]+=noiseVal;
 400               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 401               ++cs[type];
 402             }
 403           }
 404         }
 405         //whole tpc
 406         Int_t type=0;
 407         vNoiseMean[type]+=noiseVal;
 408         vNoiseRMS[type]+=noiseVal*noiseVal;
 409         ++c[type];
 410         if (!masksen){
 411           vNoiseMeanSenRegions[type]+=noiseVal;
 412           vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 413           ++cs[type];
 414         }
 415       }//end loop pads
 416     }//end loop rows
 417   }//end loop sectors (rocs)
 418
 419   //calculate mean and RMS
 420   const Double_t verySmall=0.0000000001;
 421   for (UInt_t i=0;i<infoSize;++i){
 422     Double_t mean=0;
 423     Double_t rms=0;
 424     Double_t meanSen=0;
 425     Double_t rmsSen=0;
 426
 427     if (c[i]>verySmall){
 428 //       printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
 429       mean=vNoiseMean[i]/c[i];
 430       rms=vNoiseRMS[i];
 431       rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
 432     }
 433     vNoiseMean[i]=mean;
 434     vNoiseRMS[i]=rms;
 435
 436     if (cs[i]>verySmall){
 437       meanSen=vNoiseMeanSenRegions[i]/cs[i];
 438       rmsSen=vNoiseRMSSenRegions[i];
 439       rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
 440     }
 441     vNoiseMeanSenRegions[i]=meanSen;
 442     vNoiseRMSSenRegions[i]=rmsSen;
 443   }
 444 }
 445
 446 //_____________________________________________________________________________________
 447 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
 448 {
 449   //
 450   // Process the Pulser information
 451   // vMeanTime:     pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
 452   //
 453
 454   const UInt_t infoSize=4;
 455   //reset counters to error number
 456   vMeanTime.ResizeTo(infoSize);
 457   vMeanTime.Zero();
 458   //counter
 459   TVectorD c(infoSize);
 460   //retrieve pulser and ALTRO data
 461   if (!fPulserTmean) return;
 462   //
 463   //get Outliers
 464   AliTPCCalROC *rocOut=0x0;
 465   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 466     AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
 467     if (!tmeanROC) continue;
 468     rocOut=fPulserOutlier->GetCalROC(isec);
 469     UInt_t nchannels=tmeanROC->GetNchannels();
 470     for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
 471       if (rocOut && rocOut->GetValue(ichannel)) continue;
 472       Float_t val=tmeanROC->GetValue(ichannel);
 473       Int_t type=isec/18;
 474       vMeanTime[type]+=val;
 475       ++c[type];
 476     }
 477   }
 478   //calculate mean
 479   for (UInt_t itype=0; itype<infoSize; ++itype){
 480     if (c[itype]>0) vMeanTime[itype]/=c[itype];
 481     else vMeanTime[itype]=0;
 482   }
 483 }
 484 //_____________________________________________________________________________________
 485 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
 486 {
 487   //
 488   // Get Values from ALTRO configuration data
 489   //
 490   nMasked=-1;
 491   if (!fALTROMasked) return;
 492   nMasked=0;
 493   for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
 494     AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
 495     for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
 496       if (rocMasked->GetValue(ichannel)) ++nMasked;
 497     }
 498   }
 499 }
 500 //_____________________________________________________________________________________
 501 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
 502 {
 503   //
 504   // Proces Goofie values, return statistical information of the currently set goofieArray
 505   // The meaning of the entries are given below
 506   /*
 507   1       TPC_ANODE_I_A00_STAT
 508   2       TPC_DVM_CO2
 509   3       TPC_DVM_DriftVelocity
 510   4       TPC_DVM_FCageHV
 511   5       TPC_DVM_GainFar
 512   6       TPC_DVM_GainNear
 513   7       TPC_DVM_N2
 514   8       TPC_DVM_NumberOfSparks
 515   9       TPC_DVM_PeakAreaFar
 516   10      TPC_DVM_PeakAreaNear
 517   11      TPC_DVM_PeakPosFar
 518   12      TPC_DVM_PeakPosNear
 519   13      TPC_DVM_PickupHV
 520   14      TPC_DVM_Pressure
 521   15      TPC_DVM_T1_Over_P
 522   16      TPC_DVM_T2_Over_P
 523   17      TPC_DVM_T_Over_P
 524   18      TPC_DVM_TemperatureS1
 525    */
 526   if (!fGoofieArray){
 527     Int_t nsensors=19;
 528     vecEntries.ResizeTo(nsensors);
 529     vecMedian.ResizeTo(nsensors);
 530     vecMean.ResizeTo(nsensors);
 531     vecRMS.ResizeTo(nsensors);
 532     vecEntries.Zero();
 533     vecMedian.Zero();
 534     vecMean.Zero();
 535     vecRMS.Zero();
 536     return;
 537   }
 538   Double_t kEpsilon=0.0000000001;
 539   Double_t kBig=100000000000.;
 540   Int_t nsensors = fGoofieArray->NumSensors();
 541   vecEntries.ResizeTo(nsensors);
 542   vecMedian.ResizeTo(nsensors);
 543   vecMean.ResizeTo(nsensors);
 544   vecRMS.ResizeTo(nsensors);
 545   TVectorF values;
 546   for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
 547     AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
 548     if (gsensor &&  gsensor->GetGraph()){
 549       Int_t npoints = gsensor->GetGraph()->GetN();
 550       // filter zeroes
 551       values.ResizeTo(npoints);
 552       Int_t nused =0;
 553       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 554         if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
 555             TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
 556               values[nused]=gsensor->GetGraph()->GetY()[ipoint];
 557               nused++;
 558             }
 559       }
 560       //
 561       vecEntries[isensor]= nused;
 562       if (nused>1){
 563         vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
 564         vecMean[isensor]   = TMath::Mean(nused,values.GetMatrixArray());
 565         vecRMS[isensor]    = TMath::RMS(nused,values.GetMatrixArray());
 566       }
 567     }
 568   }
 569 }
 570 //_____________________________________________________________________________________
 571 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
 572 {
 573   //
 574   // check the variations of the pedestal data to the reference pedestal data
 575   // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
 576   //
 577   const Int_t npar=4;
 578   TVectorF vThres(npar); //thresholds
 579   Int_t nActive=0;       //number of active channels
 580
 581   //reset and set thresholds
 582   pedestalDeviations.ResizeTo(npar);
 583   for (Int_t i=0;i<npar;++i){
 584     pedestalDeviations.GetMatrixArray()[i]=0;
 585     vThres.GetMatrixArray()[i]=(i+1)*.5;
 586   }
 587   //check all needed data is available
 588   if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
 589   //loop over all channels
 590   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 591     AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
 592     AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
 593     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 594     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 595     UInt_t nrows=mROC->GetNrows();
 596     for (UInt_t irow=0;irow<nrows;++irow){
 597       UInt_t npads=mROC->GetNPads(irow);
 598       for (UInt_t ipad=0;ipad<npads;++ipad){
 599         //don't use masked channels;
 600         if (mROC   ->GetValue(irow,ipad)) continue;
 601         if (mRefROC->GetValue(irow,ipad)) continue;
 602         Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
 603         for (Int_t i=0;i<npar;++i){
 604           if (deviation>vThres[i])
 605             ++pedestalDeviations.GetMatrixArray()[i];
 606         }
 607         ++nActive;
 608       }//end ipad
 609     }//ind irow
 610   }//end isec
 611   if (nActive>0){
 612     for (Int_t i=0;i<npar;++i){
 613       pedestalDeviations.GetMatrixArray()[i]/=nActive;
 614     }
 615   }
 616 }
 617 //_____________________________________________________________________________________
 618 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
 619 {
 620   //
 621   // check the variations of the noise data to the reference noise data
 622   // thresholds are 5, 10, 15 and 20 percent respectively.
 623   //
 624   const Int_t npar=4;
 625   TVectorF vThres(npar); //thresholds
 626   Int_t nActive=0;       //number of active channels
 627
 628   //reset and set thresholds
 629   noiseDeviations.ResizeTo(npar);
 630   for (Int_t i=0;i<npar;++i){
 631     noiseDeviations.GetMatrixArray()[i]=0;
 632     vThres.GetMatrixArray()[i]=(i+1)*.05;
 633   }
 634   //check all needed data is available
 635   if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
 636   //loop over all channels
 637   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 638     AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
 639     AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
 640     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 641     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 642     UInt_t nrows=mROC->GetNrows();
 643     for (UInt_t irow=0;irow<nrows;++irow){
 644       UInt_t npads=mROC->GetNPads(irow);
 645       for (UInt_t ipad=0;ipad<npads;++ipad){
 646         //don't use masked channels;
 647         if (mROC   ->GetValue(irow,ipad)) continue;
 648         if (mRefROC->GetValue(irow,ipad)) continue;
 649         Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
 650         for (Int_t i=0;i<npar;++i){
 651           if (deviation>vThres[i])
 652             ++noiseDeviations.GetMatrixArray()[i];
 653         }
 654         ++nActive;
 655       }//end ipad
 656     }//ind irow
 657   }//end isec
 658   if (nActive>0){
 659     for (Int_t i=0;i<npar;++i){
 660       noiseDeviations.GetMatrixArray()[i]/=nActive;
 661     }
 662   }
 663 }
 664 //_____________________________________________________________________________________
 665 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
 666                                                 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
 667 {
 668   //
 669   // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
 670   // thresholds are .5, 1, 5 and 10 percent respectively.
 671   //
 672   //
 673   const Int_t npar=4;
 674   TVectorF vThres(npar); //thresholds
 675   Int_t nActive=0;       //number of active channels
 676
 677   //reset and set thresholds
 678   pulserQdeviations.ResizeTo(npar);
 679   for (Int_t i=0;i<npar;++i){
 680     pulserQdeviations.GetMatrixArray()[i]=0;
 681   }
 682   npadsOutOneTB=0;
 683   npadsOffAdd=0;
 684   varQMean=0;
 685   vThres.GetMatrixArray()[0]=.005;
 686   vThres.GetMatrixArray()[1]=.01;
 687   vThres.GetMatrixArray()[2]=.05;
 688   vThres.GetMatrixArray()[3]=.1;
 689   //check all needed data is available
 690   if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
 691   //
 692   UpdateRefPulserOutlierMap();
 693   //loop over all channels
 694   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 695     AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
 696     AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
 697     AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
 698 //     AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
 699     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 700     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 701     AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
 702     Float_t pt_mean=ptROC->GetMean(oROC);
 703     UInt_t nrows=mROC->GetNrows();
 704     for (UInt_t irow=0;irow<nrows;++irow){
 705       UInt_t npads=mROC->GetNPads(irow);
 706       for (UInt_t ipad=0;ipad<npads;++ipad){
 707         //don't use masked channels;
 708         if (mROC   ->GetValue(irow,ipad)) continue;
 709         if (mRefROC->GetValue(irow,ipad)) continue;
 710         //don't user edge pads
 711         if (ipad==0||ipad==npads-1) continue;
 712         //data
 713         Float_t pq=pqROC->GetValue(irow,ipad);
 714         Float_t pqRef=pqRefROC->GetValue(irow,ipad);
 715         Float_t pt=ptROC->GetValue(irow,ipad);
 716 //         Float_t ptRef=ptRefROC->GetValue(irow,ipad);
 717         //comparisons q
 718         Float_t deviation=TMath::Abs(pq/pqRef-1);
 719         for (Int_t i=0;i<npar;++i){
 720           if (deviation>vThres[i])
 721             ++pulserQdeviations.GetMatrixArray()[i];
 722         }
 723         if (pqRef>11&&pq<11) ++npadsOffAdd;
 724         varQMean+=pq-pqRef;
 725         //comparisons t
 726         if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
 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       pulserQdeviations.GetMatrixArray()[i]/=nActive;
 734       varQMean/=nActive;
 735     }
 736   }
 737 }
 738 //_____________________________________________________________________________________
 739 void AliTPCcalibDButil::UpdatePulserOutlierMap()
 740 {
 741   //
 742   //
 743   //
 744   PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
 745 }
 746 //_____________________________________________________________________________________
 747 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
 748 {
 749   //
 750   //
 751   //
 752   PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
 753 }
 754 //_____________________________________________________________________________________
 755 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
 756 {
 757   //
 758   // Create a map that contains outliers from the Pulser calibration data.
 759   // The outliers include masked channels, edge pads and pads with
 760   //   too large timing and charge variations.
 761   // fNpulserOutliers is the number of outliers in the Pulser calibration data.
 762   //   those do not contain masked and edge pads
 763   //
 764   if (!pulT||!pulQ) {
 765     //reset map
 766     pulOut->Multiply(0.);
 767     fNpulserOutliers=-1;
 768     return;
 769   }
 770   AliTPCCalROC *rocMasked=0x0;
 771   fNpulserOutliers=0;
 772
 773   //Create Outlier Map
 774   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 775     AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
 776     AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
 777     AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
 778     if (!tmeanROC||!qmeanROC) {
 779       //reset outliers in this ROC
 780       outROC->Multiply(0.);
 781       continue;
 782     }
 783     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 784 //     Double_t dummy=0;
 785 //     Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
 786 //     Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
 787     UInt_t nrows=tmeanROC->GetNrows();
 788     for (UInt_t irow=0;irow<nrows;++irow){
 789       UInt_t npads=tmeanROC->GetNPads(irow);
 790       for (UInt_t ipad=0;ipad<npads;++ipad){
 791         Int_t outlier=0,masked=0;
 792         Float_t q=qmeanROC->GetValue(irow,ipad);
 793         Float_t t=tmeanROC->GetValue(irow,ipad);
 794         //masked channels are outliers
 795         if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
 796         //edge pads are outliers
 797         if (ipad==0||ipad==npads-1) masked=1;
 798         //channels with too large charge or timing deviation from the meadian are outliers
 799 //         if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
 800         if (q<fPulQminLimit && !masked) outlier=1;
 801         //check for nan
 802         if ( !(q<10000000) || !(t<10000000)) outlier=1;
 803         outROC->SetValue(irow,ipad,outlier+masked);
 804         fNpulserOutliers+=outlier;
 805       }
 806     }
 807   }
 808 }
 809 //_____________________________________________________________________________________
 810 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model, Double_t &gyA, Double_t &gyC, Double_t &chi2A, Double_t &chi2C )
 811 {
 812   //
 813   // Create pad time0 object from pulser and/or CE data, depending on the selected model
 814   // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
 815   // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
 816   // Model 2: use CE data and a combination CE fit + pulser in the outlier regions.
 817   //
 818   // In case model 2 is invoked - gy arival time gradient is also returned
 819   //
 820   gyA=0;
 821   gyC=0;
 822   AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
 823   // decide between different models
 824   if (model==0||model==1){
 825     TVectorD vMean;
 826     if (model==1) ProcessPulser(vMean);
 827     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 828       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 829       if (!rocPulTmean) continue;
 830       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 831       AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
 832       Float_t mean=rocPulTmean->GetMean(rocOut);
 833       //treat case where a whole partition is masked
 834       if (mean==0) mean=rocPulTmean->GetMean();
 835       if (model==1) {
 836         Int_t type=isec/18;
 837         mean=vMean[type];
 838       }
 839       UInt_t nrows=rocTime0->GetNrows();
 840       for (UInt_t irow=0;irow<nrows;++irow){
 841         UInt_t npads=rocTime0->GetNPads(irow);
 842         for (UInt_t ipad=0;ipad<npads;++ipad){
 843           Float_t time=rocPulTmean->GetValue(irow,ipad);
 844           //in case of an outlier pad use the mean of the altro values.
 845           //This should be the most precise guess in that case.
 846           if (rocOut->GetValue(irow,ipad)) {
 847             time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
 848             if (time==0) time=mean;
 849           }
 850           Float_t val=time-mean;
 851           rocTime0->SetValue(irow,ipad,val);
 852         }
 853       }
 854     }
 855   } else if (model==2){
 856     Double_t pgya,pgyc,pchi2a,pchi2c;
 857     AliTPCCalPad * padPulser = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
 858     fCETmean->Add(padPulser,-1.);
 859     TVectorD vA,vC;
 860     AliTPCCalPad outCE("outCE","outCE");
 861     Int_t nOut;
 862     ProcessCEdata("(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)++(ly/lx)^2",vA,vC,nOut,chi2A, chi2C,&outCE);
 863     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++0++gy++0++(lx-134)++0++0",vA,vC);
 864 //     AliTPCCalPad *padFit=AliTPCCalPad::CreateCalPadFit("1++(sector<36)++gy++gx++(lx-134)++(sector<36)*(lx-134)",vA,vC);
 865     if (!padFit) { delete padPulser; return 0;}
 866     gyA=vA[2];
 867     gyC=vC[2];
 868     fCETmean->Add(padPulser,1.);
 869     padTime0->Add(fCETmean);
 870     padTime0->Add(padFit,-1);
 871     delete padPulser;
 872     TVectorD vFitROC;
 873     TMatrixD mFitROC;
 874     Float_t chi2;
 875     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 876       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 877       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 878       AliTPCCalROC *rocOutPul=fPulserOutlier->GetCalROC(isec);
 879       AliTPCCalROC *rocOutCE=outCE.GetCalROC(isec);
 880       rocTime0->GlobalFit(rocOutCE,kFALSE,vFitROC,mFitROC,chi2);
 881       AliTPCCalROC *rocCEfit=AliTPCCalROC::CreateGlobalFitCalROC(vFitROC, isec);
 882       Float_t mean=rocPulTmean->GetMean(rocOutPul);
 883       if (mean==0) mean=rocPulTmean->GetMean();
 884       UInt_t nrows=rocTime0->GetNrows();
 885       for (UInt_t irow=0;irow<nrows;++irow){
 886         UInt_t npads=rocTime0->GetNPads(irow);
 887         for (UInt_t ipad=0;ipad<npads;++ipad){
 888           Float_t timePulser=rocPulTmean->GetValue(irow,ipad)-mean;
 889           if (rocOutCE->GetValue(irow,ipad)){
 890             Float_t valOut=rocCEfit->GetValue(irow,ipad);
 891             if (!rocOutPul->GetValue(irow,ipad)) valOut+=timePulser;
 892             rocTime0->SetValue(irow,ipad,valOut);
 893           }
 894         }
 895       }
 896       delete rocCEfit;
 897     }
 898     delete padFit;
 899   }
 900   Double_t median = padTime0->GetMedian();
 901   padTime0->Add(-median);  // normalize to median
 902   return padTime0;
 903 }
 904 //_____________________________________________________________________________________
 905 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
 906 {
 907   if (roc==0) return 0.;
 908   const Int_t sector=roc->GetSector();
 909   AliTPCROC *tpcRoc=AliTPCROC::Instance();
 910   const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
 911   Float_t mean=0;
 912   Int_t   n=0;
 913
 914   //loop over a small range around the requested pad (+-10 rows/pads)
 915   for (Int_t irow=row-10;irow<row+10;++irow){
 916     if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
 917     for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
 918       if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
 919       const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
 920       if (altroRoc!=altroCurr) continue;
 921       if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
 922       Float_t val=roc->GetValue(irow,ipad);
 923       mean+=val;
 924       ++n;
 925     }
 926   }
 927   if (n>0) mean/=n;
 928   return mean;
 929 }
 930 //_____________________________________________________________________________________
 931 void AliTPCcalibDButil::SetRefFile(const char* filename)
 932 {
 933   //
 934   // load cal pad objects form the reference file
 935   //
 936   TDirectory *currDir=gDirectory;
 937   TFile f(filename);
 938   fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
 939   fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
 940   //pulser data
 941   fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
 942   fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
 943   fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
 944   //CE data
 945   fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
 946   fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
 947   fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
 948   //Altro data
 949 //   fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
 950 //   fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
 951 //   fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
 952 //   fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
 953   fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
 954   f.Close();
 955   currDir->cd();
 956 }
 957
 958
 959
 960
 961 AliTPCCalPad *AliTPCcalibDButil::CreateCEOutlyerMap( Int_t & noutliersCE, AliTPCCalPad *ceOut, Float_t minSignal, Float_t cutTrmsMin,  Float_t cutTrmsMax, Float_t cutMaxDistT){
 962   //
 963   // Author:  marian.ivanov@cern.ch
 964   //
 965   // Create outlier map for CE study
 966   // Parameters:
 967   //  Return value - outlyer map
 968   //  noutlyersCE  - number of outlyers
 969   //  minSignal    - minimal total Q signal
 970   //  cutRMSMin    - minimal width of the signal in respect to the median
 971   //  cutRMSMax    - maximal width of the signal in respect to the median
 972   //  cutMaxDistT  - maximal deviation from time median per chamber
 973   //
 974   // Outlyers criteria:
 975   // 0. Exclude masked pads
 976   // 1. Exclude first two rows in IROC and last two rows in OROC
 977   // 2. Exclude edge pads
 978   // 3. Exclude channels with too large variations
 979   // 4. Exclude pads with too small signal
 980   // 5. Exclude signal with outlyers RMS
 981   // 6. Exclude channels to far from the chamber median
 982   noutliersCE=0;
 983   //create outlier map
 984   AliTPCCalPad *out=ceOut;
 985   if (!out)     out= new AliTPCCalPad("outCE","outCE");
 986   AliTPCCalROC *rocMasked=0x0;
 987   if (!fCETmean) return 0;
 988   if (!fCETrms) return 0;
 989   if (!fCEQmean) return 0;
 990   //
 991   //loop over all channels
 992   //
 993   Double_t rmsMedian         = fCETrms->GetMedian();
 994   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 995     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 996     if (fALTROMasked) rocMasked= fALTROMasked->GetCalROC(iroc);
 997     AliTPCCalROC *rocOut       = out->GetCalROC(iroc);
 998     AliTPCCalROC *rocCEQ       = fCEQmean->GetCalROC(iroc);
 999     AliTPCCalROC *rocCETrms    = fCETrms->GetCalROC(iroc);
1000     Double_t trocMedian        = rocData->GetMedian();
1001     //
1002     if (!rocData) {
1003       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
1004       rocOut->Add(1.);
1005       continue;
1006     }
1007     //
1008     //select outliers
1009     UInt_t nrows=rocData->GetNrows();
1010     for (UInt_t irow=0;irow<nrows;++irow){
1011       UInt_t npads=rocData->GetNPads(irow);
1012       for (UInt_t ipad=0;ipad<npads;++ipad){
1013         rocOut->SetValue(irow,ipad,0);
1014         Float_t valTmean=rocData->GetValue(irow,ipad);
1015         Float_t valQmean=rocCEQ->GetValue(irow,ipad);
1016         Float_t valTrms =rocCETrms->GetValue(irow,ipad);
1017         //0. exclude masked pads
1018         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
1019           rocOut->SetValue(irow,ipad,1);
1020           continue;
1021         }
1022         //1. exclude first two rows in IROC and last two rows in OROC
1023         if (iroc<36){
1024           if (irow<2) rocOut->SetValue(irow,ipad,1);
1025         } else {
1026           if (irow>nrows-3) rocOut->SetValue(irow,ipad,1);
1027         }
1028         //2. exclude edge pads
1029         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
1030         //exclude values that are exactly 0
1031         if (valTmean==0) {
1032           rocOut->SetValue(irow,ipad,1);
1033           ++noutliersCE;
1034         }
1035         //3.  exclude channels with too large variations
1036         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
1037           rocOut->SetValue(irow,ipad,1);
1038           ++noutliersCE;
1039         }
1040         //
1041         //4.  exclude channels with too small signal
1042         if (valQmean<minSignal) {
1043           rocOut->SetValue(irow,ipad,1);
1044           ++noutliersCE;
1045         }
1046         //
1047         //5. exclude channels with too small rms
1048         if (valTrms<cutTrmsMin*rmsMedian || valTrms>cutTrmsMax*rmsMedian){
1049           rocOut->SetValue(irow,ipad,1);
1050           ++noutliersCE;
1051         }
1052         //
1053         //6. exclude channels to far from the chamber median
1054         if (TMath::Abs(valTmean-trocMedian)>cutMaxDistT){
1055           rocOut->SetValue(irow,ipad,1);
1056           ++noutliersCE;
1057         }
1058       }
1059     }
1060   }
1061   //
1062   return out;
1063 }
1064
1065
1066 AliTPCCalPad *AliTPCcalibDButil::CreatePulserOutlyerMap(Int_t &noutliersPulser, AliTPCCalPad *pulserOut,Float_t cutTime, Float_t cutnRMSQ, Float_t cutnRMSrms){
1067   //
1068   // Author: marian.ivanov@cern.ch
1069   //
1070   // Create outlier map for Pulser
1071   // Parameters:
1072   //  Return value     - outlyer map
1073   //  noutlyersPulser  - number of outlyers
1074   //  cutTime          - absolute cut - distance to the median of chamber
1075   //  cutnRMSQ         - nsigma cut from median  q distribution per chamber
1076   //  cutnRMSrms       - nsigma cut from median  rms distribution
1077   // Outlyers criteria:
1078   // 0. Exclude masked pads
1079   // 1. Exclude time outlyers (default 3 time bins)
1080   // 2. Exclude q outlyers    (default 5 sigma)
1081   // 3. Exclude rms outlyers  (default 5 sigma)
1082   noutliersPulser=0;
1083   AliTPCCalPad *out=pulserOut;
1084   if (!out)     out= new AliTPCCalPad("outPulser","outPulser");
1085   AliTPCCalROC *rocMasked=0x0;
1086   if (!fPulserTmean) return 0;
1087   if (!fPulserTrms) return 0;
1088   if (!fPulserQmean) return 0;
1089   //
1090   //loop over all channels
1091   //
1092   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1093     if (fALTROMasked)   rocMasked= fALTROMasked->GetCalROC(iroc);
1094     AliTPCCalROC *rocData       = fPulserTmean->GetCalROC(iroc);
1095     AliTPCCalROC *rocOut        = out->GetCalROC(iroc);
1096     AliTPCCalROC *rocPulserQ    = fPulserQmean->GetCalROC(iroc);
1097     AliTPCCalROC *rocPulserTrms = fPulserTrms->GetCalROC(iroc);
1098     //
1099     Double_t rocMedianT         = rocData->GetMedian();
1100     Double_t rocMedianQ         = rocPulserQ->GetMedian();
1101     Double_t rocRMSQ            = rocPulserQ->GetRMS();
1102     Double_t rocMedianTrms      = rocPulserTrms->GetMedian();
1103     Double_t rocRMSTrms         = rocPulserTrms->GetRMS();
1104     for (UInt_t ichannel=0;ichannel<rocData->GetNchannels();++ichannel){
1105       rocOut->SetValue(ichannel,0);
1106       Float_t valTmean=rocData->GetValue(ichannel);
1107       Float_t valQmean=rocPulserQ->GetValue(ichannel);
1108       Float_t valTrms =rocPulserTrms->GetValue(ichannel);
1109       Int_t isOut=0;
1110       if (TMath::Abs(valTmean-rocMedianT)>cutTime) isOut=1;
1111       if (TMath::Abs(valQmean-rocMedianQ)>cutnRMSQ*rocRMSQ) isOut=1;
1112       if (TMath::Abs(valTrms-rocMedianTrms)>cutnRMSrms*rocRMSTrms) isOut=1;
1113       rocOut->SetValue(ichannel,isOut);
1114       if (isOut) noutliersPulser++;
1115     }
1116   }
1117   return out;
1118 }
1119
1120
1121 AliTPCCalPad *AliTPCcalibDButil::CreatePadTime0CE(TVectorD &fitResultsA, TVectorD&fitResultsC, Int_t &nOut, Double_t &chi2A, Double_t &chi2C, const char *dumpfile){
1122   //
1123   // Author : Marian Ivanov
1124   // Create pad time0 correction map using information from the CE and from pulser
1125   //
1126   //
1127   // Return PadTime0 to be used for time0 relative alignment
1128   // if dump file specified intermediat results are dumped to the fiel and can be visualized
1129   // using $ALICE_ROOT/TPC/script/gui application
1130   //
1131   // fitResultsA - fitParameters A side
1132   // fitResultsC - fitParameters C side
1133   // chi2A       - chi2/ndf for A side (assuming error 1 time bin)
1134   // chi2C       - chi2/ndf for C side (assuming error 1 time bin)
1135   //
1136   //
1137   // Algorithm:
1138   // 1. Find outlier map for CE
1139   // 2. Find outlier map for Pulser
1140   // 3. Replace outlier by median at given sector  (median without outliers)
1141   // 4. Substract from the CE data pulser
1142   // 5. Fit the CE with formula
1143   //    5.1) (IROC-OROC) offset
1144   //    5.2) gx
1145   //    5.3) gy
1146   //    5.4) (lx-xmid)
1147   //    5.5) (IROC-OROC)*(lx-xmid)
1148   //    5.6) (ly/lx)^2
1149   // 6. Substract gy fit dependence from the CE data
1150   // 7. Add pulser back to CE data
1151   // 8. Replace outliers by fit value - median of diff per given chamber -GY fit
1152   // 9. return CE data
1153   //
1154   // Time0 <= padCE = padCEin  -padCEfitGy  - if not outlier
1155   // Time0 <= padCE = padFitAll-padCEfitGy  - if outlier
1156
1157   // fit formula
1158   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)";
1159   // output for fit formula
1160   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)";
1161   // gy part of formula
1162   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)";
1163   //
1164   //
1165   if (!fCETmean) return 0;
1166   Double_t pgya,pgyc,pchi2a,pchi2c;
1167   AliTPCCalPad * padPulserOut = CreatePulserOutlyerMap(nOut);
1168   AliTPCCalPad * padCEOut     = CreateCEOutlyerMap(nOut);
1169
1170   AliTPCCalPad * padPulser    = CreatePadTime0(1,pgya,pgyc,pchi2a,pchi2c);
1171   AliTPCCalPad * padCE        = new AliTPCCalPad(*fCETmean);
1172   AliTPCCalPad * padCEIn      = new AliTPCCalPad(*fCETmean);
1173   AliTPCCalPad * padOut       = new AliTPCCalPad("padOut","padOut");
1174   padPulser->SetName("padPulser");
1175   padPulserOut->SetName("padPulserOut");
1176   padCE->SetName("padCE");
1177   padCEIn->SetName("padCEIn");
1178   padCEOut->SetName("padCEOut");
1179   padOut->SetName("padOut");
1180
1181   //
1182   // make combined outlyers map
1183   // and replace outlyers in maps with median for chamber
1184   //
1185   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1186     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1187     AliTPCCalROC * rocPulser    = padPulser->GetCalROC(iroc);
1188     AliTPCCalROC * rocPulserOut = padPulserOut->GetCalROC(iroc);
1189     AliTPCCalROC * rocCEOut     = padCEOut->GetCalROC(iroc);
1190     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1191     Double_t ceMedian           = rocCE->GetMedian(rocCEOut);
1192     Double_t pulserMedian       = rocPulser->GetMedian(rocCEOut);
1193     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1194       if (rocPulserOut->GetValue(ichannel)>0) {
1195         rocPulser->SetValue(ichannel,pulserMedian);
1196         rocOut->SetValue(ichannel,1);
1197       }
1198       if (rocCEOut->GetValue(ichannel)>0) {
1199         rocCE->SetValue(ichannel,ceMedian);
1200         rocOut->SetValue(ichannel,1);
1201       }
1202     }
1203   }
1204   //
1205   // remove pulser time 0
1206   //
1207   padCE->Add(padPulser,-1);
1208   //
1209   // Make fits
1210   //
1211   TMatrixD dummy;
1212   Float_t chi2Af,chi2Cf;
1213   padCE->GlobalSidesFit(padOut,formulaIn,fitResultsA,fitResultsC,dummy,dummy,chi2Af,chi2Cf);
1214   chi2A=chi2Af;
1215   chi2C=chi2Cf;
1216   //
1217   AliTPCCalPad *padCEFitGY=AliTPCCalPad::CreateCalPadFit(formulaOut,fitResultsA,fitResultsC);
1218   padCEFitGY->SetName("padCEFitGy");
1219   //
1220   AliTPCCalPad *padCEFit  =AliTPCCalPad::CreateCalPadFit(formulaAll,fitResultsA,fitResultsC);
1221   padCEFit->SetName("padCEFit");
1222   //
1223   AliTPCCalPad* padCEDiff  = new AliTPCCalPad(*padCE);
1224   padCEDiff->SetName("padCEDiff");
1225   padCEDiff->Add(padCEFit,-1.);
1226   //
1227   //
1228   padCE->Add(padCEFitGY,-1.);
1229
1230   padCE->Add(padPulser,1.);
1231   Double_t padmedian = padCE->GetMedian();
1232   padCE->Add(-padmedian);  // normalize to median
1233   //
1234   // Replace outliers by fit value - median of diff per given chamber -GY fit
1235   //
1236   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
1237     AliTPCCalROC * rocOut       = padOut->GetCalROC(iroc);
1238     AliTPCCalROC * rocCE        = padCE->GetCalROC(iroc);
1239     AliTPCCalROC * rocCEFit     = padCEFit->GetCalROC(iroc);
1240     AliTPCCalROC * rocCEFitGY   = padCEFitGY->GetCalROC(iroc);
1241     AliTPCCalROC * rocCEDiff    = padCEDiff->GetCalROC(iroc);
1242     //
1243     Double_t diffMedian         = rocCEDiff->GetMedian(rocOut);
1244     for (UInt_t ichannel=0;ichannel<rocOut->GetNchannels();++ichannel){
1245       if (rocOut->GetValue(ichannel)==0) continue;
1246       Float_t value=rocCEFit->GetValue(ichannel)-rocCEFitGY->GetValue(ichannel)-diffMedian-padmedian;
1247       rocCE->SetValue(ichannel,value);
1248     }
1249   }
1250   //
1251   //
1252   if (dumpfile){
1253     //dump to the file - result can be visualized
1254     AliTPCPreprocessorOnline preprocesor;
1255     preprocesor.AddComponent(new AliTPCCalPad(*padCE));
1256     preprocesor.AddComponent(new AliTPCCalPad(*padCEIn));
1257     preprocesor.AddComponent(new AliTPCCalPad(*padCEFit));
1258     preprocesor.AddComponent(new AliTPCCalPad(*padOut));
1259     //
1260     preprocesor.AddComponent(new AliTPCCalPad(*padCEFitGY));
1261     preprocesor.AddComponent(new AliTPCCalPad(*padCEDiff));
1262     //
1263     preprocesor.AddComponent(new AliTPCCalPad(*padCEOut));
1264     preprocesor.AddComponent(new AliTPCCalPad(*padPulser));
1265     preprocesor.AddComponent(new AliTPCCalPad(*padPulserOut));
1266     preprocesor.DumpToFile(dumpfile);
1267   }
1268   delete padPulser;
1269   delete padPulserOut;
1270   delete padCEIn;
1271   delete padCEOut;
1272   delete padOut;
1273   delete padCEDiff;
1274   delete padCEFitGY;
1275   return padCE;
1276 }
1277
1278
1279
1280
1281
1282 Int_t AliTPCcalibDButil::GetNearest(TGraph *graph, Double_t xref, Double_t &dx, Double_t &y){
1283   //
1284   // find the closest point to xref  in x  direction
1285   // return dx and value
1286   Int_t index=0;
1287   index = TMath::BinarySearch(graph->GetN(), graph->GetX(),xref);
1288   if (index<0) index=0;
1289   if (index>=graph->GetN()-1) index=graph->GetN()-2;
1290   if (xref-graph->GetX()[index]>graph->GetX()[index]-xref) index++;
1291   dx = xref-graph->GetX()[index];
1292   y  = graph->GetY()[index];
1293   return index;
1294 }
1295
1296
1297 Double_t  AliTPCcalibDButil::GetTriggerOffsetTPC(Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1298   //
1299   // Get the correction of the trigger offset
1300   // combining information from the laser track calibration
1301   // and from cosmic calibration
1302   //
1303   // run       - run number
1304   // timeStamp - tim stamp in seconds
1305   // deltaT    - integration period to calculate offset
1306   // deltaTLaser -max validity of laser data
1307   // valType   - 0 - median, 1- mean
1308   //
1309   // Integration vaues are just recomendation - if not possible to get points
1310   // automatically increase the validity by factor 2
1311   // (recursive algorithm until one month of data taking)
1312   //
1313   //
1314   const Float_t kLaserCut=0.0005;
1315   const Int_t   kMaxPeriod=3600*24*30*3; // 3 month max
1316   const Int_t   kMinPoints=20;
1317   //
1318   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1319   if (!array) {
1320     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1321   }
1322   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1323   if (!array) return 0;
1324   //
1325   TGraphErrors *laserA[3]={0,0,0};
1326   TGraphErrors *laserC[3]={0,0,0};
1327   TGraphErrors *cosmicAll=0;
1328   laserA[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1329   laserC[1]=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1330   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1331   //
1332   //
1333   if (!cosmicAll) return 0;
1334   Int_t nmeasC=cosmicAll->GetN();
1335   Float_t *tdelta = new Float_t[nmeasC];
1336   Int_t nused=0;
1337   for (Int_t i=0;i<nmeasC;i++){
1338     if (TMath::Abs(cosmicAll->GetX()[i]-timeStamp)>deltaT) continue;
1339     Float_t ccosmic=cosmicAll->GetY()[i];
1340     Double_t yA=0,yC=0,dA=0,dC=0;
1341     if (laserA[1]) GetNearest(laserA[1], cosmicAll->GetX()[i],dA,yA);
1342     if (laserC[1]) GetNearest(laserC[1], cosmicAll->GetX()[i],dC,yC);
1343     //yA=laserA[1]->Eval(cosmicAll->GetX()[i]);
1344     //yC=laserC[1]->Eval(cosmicAll->GetX()[i]);
1345     //
1346     if (TMath::Sqrt(dA*dA+dC*dC)>deltaTLaser) continue;
1347     Float_t claser=0;
1348     if (TMath::Abs(yA-yC)<kLaserCut) {
1349       claser=(yA-yC)*0.5;
1350     }else{
1351       if (i%2==0)  claser=yA;
1352       if (i%2==1)  claser=yC;
1353     }
1354     tdelta[nused]=ccosmic-claser;
1355     nused++;
1356   }
1357   if (nused<kMinPoints &&deltaT<kMaxPeriod) return  AliTPCcalibDButil::GetTriggerOffsetTPC(run, timeStamp, deltaT*2,deltaTLaser);
1358   Double_t median = TMath::Median(nused,tdelta);
1359   Double_t mean  = TMath::Mean(nused,tdelta);
1360   delete tdelta;
1361   return (valType==0) ? median:mean;
1362 }
1363
1364 Double_t  AliTPCcalibDButil::GetVDriftTPC(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Double_t deltaTLaser, Int_t valType){
1365   //
1366   // Get the correction of the drift velocity
1367   // combining information from the laser track calibration
1368   // and from cosmic calibration
1369   //
1370   // dist      - return value - distance to closest point in graph
1371   // run       - run number
1372   // timeStamp - tim stamp in seconds
1373   // deltaT    - integration period to calculate time0 offset
1374   // deltaTLaser -max validity of laser data
1375   // valType   - 0 - median, 1- mean
1376   //
1377   // Integration vaues are just recomendation - if not possible to get points
1378   // automatically increase the validity by factor 2
1379   // (recursive algorithm until one month of data taking)
1380   //
1381   //
1382   //
1383   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1384   if (!array) {
1385     AliTPCcalibDB::Instance()->UpdateRunInformations(run,kFALSE);
1386   }
1387   array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1388   if (!array) return 0;
1389   TGraphErrors *cosmicAll=0;
1390   cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1391   if (!cosmicAll) return 0;
1392   Double_t grY=0;
1393   AliTPCcalibDButil::GetNearest(cosmicAll,timeStamp,dist,grY);
1394
1395   Double_t t0= AliTPCcalibDButil::GetTriggerOffsetTPC(run,timeStamp, deltaT, deltaTLaser,valType);
1396   Double_t vcosmic=  AliTPCcalibDButil::EvalGraphConst(cosmicAll, timeStamp);
1397   if (timeStamp>cosmicAll->GetX()[cosmicAll->GetN()-1])  vcosmic=cosmicAll->GetY()[cosmicAll->GetN()-1];
1398   if (timeStamp<cosmicAll->GetX()[0])  vcosmic=cosmicAll->GetY()[0];
1399   return  vcosmic+t0;
1400
1401   /*
1402     Example usage:
1403
1404     Int_t run=89000
1405     TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1406     cosmicAll =(TGraphErrors*)array->FindObject("TGRAPHERRORS_MEAN_VDRIFT_COSMICS_ALL");
1407     laserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1408     //
1409     Double_t *yvd= new Double_t[cosmicAll->GetN()];
1410     Double_t *yt0= new Double_t[cosmicAll->GetN()];
1411     for (Int_t i=0; i<cosmicAll->GetN();i++) yvd[i]=AliTPCcalibDButil::GetVDriftTPC(run,cosmicAll->GetX()[i]);
1412     for (Int_t i=0; i<cosmicAll->GetN();i++) yt0[i]=AliTPCcalibDButil::GetTriggerOffsetTPC(run,cosmicAll->GetX()[i]);
1413
1414     TGraph *pcosmicVd=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yvd);
1415     TGraph *pcosmicT0=new TGraph(cosmicAll->GetN(), cosmicAll->GetX(), yt0);
1416
1417   */
1418
1419 }
1420
1421 Double_t  AliTPCcalibDButil::GetVDriftTPCLaserTracks(Double_t &dist, Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1422   //
1423   // Get the correction of the drift velocity using the laser tracks calbration
1424   //
1425   // run       - run number
1426   // timeStamp - tim stamp in seconds
1427   // deltaT    - integration period to calculate time0 offset
1428   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1429   // Note in case no data form both A and C side - the value from active side used
1430   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
1431   TGraphErrors *grlaserA=0;
1432   TGraphErrors *grlaserC=0;
1433   Double_t vlaserA=0, vlaserC=0;
1434   if (!array) return 0;
1435   grlaserA=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_A");
1436   grlaserC=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DRIFT_LASER_ALL_C");
1437   Double_t deltaY;
1438   if (grlaserA) {
1439     AliTPCcalibDButil::GetNearest(grlaserA,timeStamp,dist,deltaY);
1440     if (TMath::Abs(dist)>deltaT)  vlaserA= deltaY;
1441     else  vlaserA = AliTPCcalibDButil::EvalGraphConst(grlaserA,timeStamp);
1442   }
1443   if (grlaserC) {
1444     AliTPCcalibDButil::GetNearest(grlaserC,timeStamp,dist,deltaY);
1445     if (TMath::Abs(dist)>deltaT)  vlaserC= deltaY;
1446     else  vlaserC = AliTPCcalibDButil::EvalGraphConst(grlaserC,timeStamp);
1447   }
1448   if (side==0) return vlaserA;
1449   if (side==1) return vlaserC;
1450   Double_t mdrift=(vlaserA+vlaserC)*0.5;
1451   if (!grlaserA) return vlaserC;
1452   if (!grlaserC) return vlaserA;
1453   return mdrift;
1454 }
1455
1456
1457 Double_t  AliTPCcalibDButil::GetVDriftTPCCE(Double_t &dist,Int_t run, Int_t timeStamp, Double_t deltaT, Int_t side){
1458   //
1459   // Get the correction of the drift velocity using the CE laser data
1460   // combining information from the CE,  laser track calibration
1461   // and P/T calibration
1462   //
1463   // run       - run number
1464   // timeStamp - tim stamp in seconds
1465   // deltaT    - integration period to calculate time0 offset
1466   // side      - 0 - A side,  1 - C side, 2 - mean from both sides
1467   TObjArray *arrT     =AliTPCcalibDB::Instance()->GetCErocTtime();
1468   if (!arrT) return 0;
1469   AliTPCParam *param  =AliTPCcalibDB::Instance()->GetParameters();
1470   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1471   AliTPCCalibVdrift * driftCalib = (AliTPCCalibVdrift *)cearray->FindObject("driftPTCE");
1472   //
1473   //
1474   Double_t corrPTA = 0, corrPTC=0;
1475   Double_t ltime0A = 0, ltime0C=0;
1476   Double_t gry=0;
1477   Double_t corrA=0, corrC=0;
1478   Double_t timeA=0, timeC=0;
1479   TGraph *graphA = (TGraph*)arrT->At(72);
1480   TGraph *graphC = (TGraph*)arrT->At(73);
1481   if (!graphA && !graphC) return 0.;
1482   if (graphA &&graphA->GetN()>0) {
1483     AliTPCcalibDButil::GetNearest(graphA,timeStamp,dist,gry);
1484     timeA   = AliTPCcalibDButil::EvalGraphConst(graphA,timeStamp);
1485     Int_t mtime   =TMath::Nint((graphA->GetX()[0]+graphA->GetX()[graphA->GetN()-1])*0.5);
1486     ltime0A       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1487     if (driftCalib) corrPTA =  driftCalib->GetPTRelative(timeStamp,0);
1488     corrA = (param->GetZLength(36)/(timeA*param->GetTSample()*(1.-ltime0A)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1489     corrA-=corrPTA;
1490   }
1491   if (graphC&&graphC->GetN()>0){
1492     AliTPCcalibDButil::GetNearest(graphC,timeStamp,dist,gry);
1493     timeC=AliTPCcalibDButil::EvalGraphConst(graphC,timeStamp);
1494     Int_t mtime=TMath::Nint((graphC->GetX()[0]+graphC->GetX()[graphC->GetN()-1])*0.5);
1495     ltime0C       = GetLaserTime0(run,mtime,TMath::Nint(deltaT),0);
1496     if (driftCalib) corrPTC =  driftCalib->GetPTRelative(timeStamp,0);
1497     corrC = (param->GetZLength(54)/(timeC*param->GetTSample()*(1.-ltime0C)-param->GetL1Delay()-0*param->GetZSigma()/param->GetDriftV()))/param->GetDriftV()-1;
1498     corrC-=corrPTC;
1499   }
1500
1501   if (side ==0 ) return corrA;
1502   if (side ==1 ) return corrC;
1503   Double_t corrM= (corrA+corrC)*0.5;
1504   if (!graphA) corrM=corrC;
1505   if (!graphC) corrM=corrA;
1506   return corrM;
1507 }
1508
1509
1510
1511
1512 Int_t  AliTPCcalibDButil::MakeRunList(Int_t startRun, Int_t stopRun){
1513   //
1514   // VERY obscure method - we need something in framework
1515   // Find the TPC runs with temperature OCDB entry
1516   // cache the start and end of the run
1517   //
1518   AliCDBStorage* storage = AliCDBManager::Instance()->GetSpecificStorage("TPC/Calib/Temperature");
1519   if (!storage) storage = AliCDBManager::Instance()->GetDefaultStorage();
1520   if (!storage) return 0;
1521   TString path=storage->GetURI();
1522   TString runsT;
1523   {
1524     TString command;
1525     if (path.Contains("local")){  // find the list if local system
1526       path.ReplaceAll("local://","");
1527       path+="TPC/Calib/Temperature";
1528       command=Form("ls %s  |  sed s/_/\\ /g | awk '{print \"r\"$2}'  ",path.Data());
1529     }
1530     runsT=gSystem->GetFromPipe(command);
1531   }
1532   TObjArray *arr= runsT.Tokenize("r");
1533   if (!arr) return 0;
1534   //
1535   TArrayI indexes(arr->GetEntries());
1536   TArrayI runs(arr->GetEntries());
1537   Int_t naccept=0;
1538   {for (Int_t irun=0;irun<arr->GetEntries();irun++){
1539       Int_t irunN = atoi(arr->At(irun)->GetName());
1540       if (irunN<startRun) continue;
1541       if (irunN>stopRun) continue;
1542       runs[naccept]=irunN;
1543       naccept++;
1544     }}
1545   fRuns.Set(naccept);
1546   fRunsStart.Set(fRuns.fN);
1547   fRunsStop.Set(fRuns.fN);
1548   TMath::Sort(fRuns.fN, runs.fArray, indexes.fArray,kFALSE);
1549   for (Int_t irun=0; irun<fRuns.fN; irun++)  fRuns[irun]=runs[indexes[irun]];
1550
1551   //
1552   AliCDBEntry * entry = 0;
1553   {for (Int_t irun=0;irun<fRuns.fN; irun++){
1554       entry = AliCDBManager::Instance()->Get("TPC/Calib/Temperature",fRuns[irun]);
1555       if (!entry) continue;
1556       AliTPCSensorTempArray *  tmpRun = dynamic_cast<AliTPCSensorTempArray*>(entry->GetObject());
1557       if (!tmpRun) continue;
1558       fRunsStart[irun]=tmpRun->GetStartTime().GetSec();
1559       fRunsStop[irun]=tmpRun->GetEndTime().GetSec();
1560       //      printf("irun\t%d\tRun\t%d\t%d\t%d\n",irun,fRuns[irun],tmpRun->GetStartTime().GetSec(),tmpRun->GetEndTime().GetSec());
1561     }}
1562   return fRuns.fN;
1563 }
1564
1565
1566 Int_t AliTPCcalibDButil::FindRunTPC(Int_t    itime, Bool_t debug){
1567   //
1568   // binary search - find the run for given time stamp
1569   //
1570   Int_t index0  = TMath::BinarySearch(fRuns.fN, fRunsStop.fArray,itime);
1571   Int_t index1  = TMath::BinarySearch(fRuns.fN, fRunsStart.fArray,itime);
1572   Int_t cindex  = -1;
1573   for (Int_t index=index0; index<=index1; index++){
1574     if (fRunsStart[index]<=itime && fRunsStop[index]>=itime) cindex=index;
1575     if (debug) {
1576       printf("%d\t%d\t%d\n",fRuns[index], fRunsStart[index]-itime, fRunsStop[index]-itime);
1577     }
1578   }
1579   if (cindex<0) cindex =(index0+index1)/2;
1580   if (cindex<0) {
1581     return 0;
1582   }
1583   return fRuns[cindex];
1584 }
1585
1586
1587
1588
1589
1590 TGraph* AliTPCcalibDButil::FilterGraphMedian(TGraph * graph, Float_t sigmaCut,Double_t &medianY){
1591   //
1592   // filter outlyer measurement
1593   // Only points around median +- sigmaCut filtered
1594   if (!graph) return  0;
1595   Int_t kMinPoints=2;
1596   Int_t npoints0 = graph->GetN();
1597   Int_t npoints=0;
1598   Float_t  rmsY=0;
1599   Double_t *outx=new Double_t[npoints0];
1600   Double_t *outy=new Double_t[npoints0];
1601   //
1602   //
1603   if (npoints0<kMinPoints) return 0;
1604   for (Int_t iter=0; iter<3; iter++){
1605     npoints=0;
1606     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1607       if (graph->GetY()[ipoint]==0) continue;
1608       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>sigmaCut*rmsY) continue;
1609       outx[npoints]  = graph->GetX()[ipoint];
1610       outy[npoints]  = graph->GetY()[ipoint];
1611       npoints++;
1612     }
1613     if (npoints<=1) break;
1614     medianY  =TMath::Median(npoints,outy);
1615     rmsY   =TMath::RMS(npoints,outy);
1616   }
1617   TGraph *graphOut=0;
1618   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1619   return graphOut;
1620 }
1621
1622
1623 TGraph* AliTPCcalibDButil::FilterGraphMedianAbs(TGraph * graph, Float_t cut,Double_t &medianY){
1624   //
1625   // filter outlyer measurement
1626   // Only points around median +- cut filtered
1627   if (!graph) return  0;
1628   Int_t kMinPoints=2;
1629   Int_t npoints0 = graph->GetN();
1630   Int_t npoints=0;
1631   Float_t  rmsY=0;
1632   Double_t *outx=new Double_t[npoints0];
1633   Double_t *outy=new Double_t[npoints0];
1634   //
1635   //
1636   if (npoints0<kMinPoints) return 0;
1637   for (Int_t iter=0; iter<3; iter++){
1638     npoints=0;
1639     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1640       if (graph->GetY()[ipoint]==0) continue;
1641       if (iter>0 &&TMath::Abs(graph->GetY()[ipoint]-medianY)>cut) continue;
1642       outx[npoints]  = graph->GetX()[ipoint];
1643       outy[npoints]  = graph->GetY()[ipoint];
1644       npoints++;
1645     }
1646     if (npoints<=1) break;
1647     medianY  =TMath::Median(npoints,outy);
1648     rmsY   =TMath::RMS(npoints,outy);
1649   }
1650   TGraph *graphOut=0;
1651   if (npoints>1) graphOut= new TGraph(npoints,outx,outy);
1652   return graphOut;
1653 }
1654
1655
1656
1657 TGraphErrors* AliTPCcalibDButil::FilterGraphMedianErr(TGraphErrors * graph, Float_t sigmaCut,Double_t &medianY){
1658   //
1659   // filter outlyer measurement
1660   // Only points with normalized errors median +- sigmaCut filtered
1661   //
1662   Int_t kMinPoints=10;
1663   Int_t npoints0 = graph->GetN();
1664   Int_t npoints=0;
1665   Float_t  medianErr=0, rmsErr=0;
1666   Double_t *outx=new Double_t[npoints0];
1667   Double_t *outy=new Double_t[npoints0];
1668   Double_t *erry=new Double_t[npoints0];
1669   Double_t *nerry=new Double_t[npoints0];
1670   Double_t *errx=new Double_t[npoints0];
1671   //
1672   //
1673   if (npoints0<kMinPoints) return 0;
1674   for (Int_t iter=0; iter<3; iter++){
1675     npoints=0;
1676     for (Int_t ipoint=0; ipoint<npoints0; ipoint++){
1677       nerry[npoints]  = graph->GetErrorY(ipoint);
1678       if (iter>0 &&TMath::Abs(nerry[npoints]-medianErr)>sigmaCut*rmsErr) continue;
1679       erry[npoints]  = graph->GetErrorY(ipoint);
1680       outx[npoints]  = graph->GetX()[ipoint];
1681       outy[npoints]  = graph->GetY()[ipoint];
1682       errx[npoints]  = graph->GetErrorY(ipoint);
1683       npoints++;
1684     }
1685     if (npoints==0) break;
1686     medianErr=TMath::Median(npoints,erry);
1687     medianY  =TMath::Median(npoints,outy);
1688     rmsErr   =TMath::RMS(npoints,erry);
1689   }
1690   TGraphErrors *graphOut=0;
1691   if (npoints>1) graphOut= new TGraphErrors(npoints,outx,outy,errx,erry);
1692   delete []outx;
1693   delete []outy;
1694   delete []errx;
1695   delete []erry;
1696   return graphOut;
1697 }
1698
1699
1700 void AliTPCcalibDButil::Sort(TGraph *graph){
1701   //
1702   // sort array - neccessay for approx
1703   //
1704   Int_t npoints = graph->GetN();
1705   Int_t *indexes=new Int_t[npoints];
1706   Double_t *outx=new Double_t[npoints];
1707   Double_t *outy=new Double_t[npoints];
1708   TMath::Sort(npoints, graph->GetX(),indexes,kFALSE);
1709   for (Int_t i=0;i<npoints;i++) outx[i]=graph->GetX()[indexes[i]];
1710   for (Int_t i=0;i<npoints;i++) outy[i]=graph->GetY()[indexes[i]];
1711   for (Int_t i=0;i<npoints;i++) graph->GetX()[i]=outx[i];
1712   for (Int_t i=0;i<npoints;i++) graph->GetY()[i]=outy[i];
1713
1714 }
1715 void AliTPCcalibDButil::SmoothGraph(TGraph *graph, Double_t delta){
1716   //
1717   // smmoth graph - mean on the interval
1718   //
1719   Sort(graph);
1720   Int_t npoints = graph->GetN();
1721   Double_t *outy=new Double_t[npoints];
1722
1723   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1724     Double_t lx=graph->GetX()[ipoint];
1725     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1726     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1727     if (index0<0) index0=0;
1728     if (index1>=npoints-1) index1=npoints-1;
1729     if ((index1-index0)>1){
1730       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1731     }else{
1732       outy[ipoint]=graph->GetY()[ipoint];
1733     }
1734   }
1735  //  TLinearFitter  fitter(3,"pol2");
1736 //   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1737 //     Double_t lx=graph->GetX()[ipoint];
1738 //     Int_t index0=TMath::BinarySearch(npoints, graph->GetX(),lx-delta);
1739 //     Int_t index1=TMath::BinarySearch(npoints, graph->GetX(),lx+delta);
1740 //     if (index0<0) index0=0;
1741 //     if (index1>=npoints-1) index1=npoints-1;
1742 //     fitter.ClearPoints();
1743 //     for (Int_t jpoint=0;jpoint<index1-index0; jpoint++)
1744 //     if ((index1-index0)>1){
1745 //       outy[ipoint]  = TMath::Mean(index1-index0, &(graph->GetY()[index0]));
1746 //     }else{
1747 //       outy[ipoint]=graph->GetY()[ipoint];
1748 //     }
1749 //   }
1750
1751
1752
1753   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
1754     graph->GetY()[ipoint] = outy[ipoint];
1755   }
1756   delete[] outy;
1757 }
1758
1759 Double_t AliTPCcalibDButil::EvalGraphConst(TGraph *graph, Double_t xref){
1760   //
1761   // Use constant interpolation outside of range
1762   //
1763   if (!graph) {
1764     printf("AliTPCcalibDButil::EvalGraphConst: 0 pointer\n");
1765     return 0;
1766   }
1767   if (graph->GetN()<1){
1768     printf("AliTPCcalibDButil::EvalGraphConst: Empty graph");
1769     return 0;
1770   }
1771   if (xref<graph->GetX()[0]) return graph->GetY()[0];
1772   if (xref>graph->GetX()[graph->GetN()-1]) return graph->GetY()[graph->GetN()-1];
1773   return graph->Eval( xref);
1774 }
1775
1776 Float_t AliTPCcalibDButil::FilterSensor(AliDCSSensor * sensor, Double_t ymin, Double_t ymax, Double_t maxdy,  Double_t sigmaCut){
1777   //
1778   // Filter DCS sensor information
1779   //   ymin     - minimal value
1780   //   ymax     - max value
1781   //   maxdy    - maximal deirivative
1782   //   sigmaCut - cut on values and derivative in terms of RMS distribution
1783   // Return value - accepted fraction
1784   //
1785   // Algorithm:
1786   //
1787   // 0. Calculate median and rms of values in specified range
1788   // 1. Filter out outliers - median+-sigmaCut*rms
1789   //    values replaced by median
1790   //
1791   AliSplineFit * fit    = sensor->GetFit();
1792   if (!fit) return 0.;
1793   Int_t          nknots = fit->GetKnots();
1794   if (nknots==0) {
1795     delete fit;
1796     sensor->SetFit(0);
1797     return 0;
1798   }
1799   //
1800   Double_t *yin0  = new Double_t[nknots];
1801   Double_t *yin1  = new Double_t[nknots];
1802   Int_t naccept=0;
1803
1804   for (Int_t iknot=0; iknot< nknots; iknot++){
1805     if (fit->GetY0()[iknot]>ymin && fit->GetY0()[iknot]<ymax){
1806       yin0[naccept]  = fit->GetY0()[iknot];
1807       yin1[naccept]  = fit->GetY1()[iknot];
1808       if (TMath::Abs(fit->GetY1()[iknot])>maxdy) yin1[naccept]=0;
1809       naccept++;
1810     }
1811   }
1812   if (naccept<1) {
1813     delete fit;
1814     sensor->SetFit(0);
1815     return 0.;
1816   }
1817
1818   Double_t medianY0=0, medianY1=0;
1819   Double_t rmsY0   =0, rmsY1=0;
1820   medianY0 = TMath::Median(naccept, yin0);
1821   medianY1 = TMath::Median(naccept, yin1);
1822   rmsY0    = TMath::RMS(naccept, yin0);
1823   rmsY1    = TMath::RMS(naccept, yin1);
1824   naccept=0;
1825   //
1826   // 1. Filter out outliers - median+-sigmaCut*rms
1827   //    values replaced by median
1828   //    if replaced the derivative set to 0
1829   //
1830   for (Int_t iknot=0; iknot< nknots; iknot++){
1831     Bool_t isOK=kTRUE;
1832     if (TMath::Abs(fit->GetY0()[iknot]-medianY0)>sigmaCut*rmsY0) isOK=kFALSE;
1833     if (TMath::Abs(fit->GetY1()[iknot]-medianY1)>sigmaCut*rmsY1) isOK=kFALSE;
1834     if (nknots<2) fit->GetY1()[iknot]=0;
1835     if (TMath::Abs(fit->GetY1()[iknot])>maxdy) fit->GetY1()[iknot]=0;
1836     if (!isOK){
1837       fit->GetY0()[iknot]=medianY0;
1838       fit->GetY1()[iknot]=0;
1839     }else{
1840       naccept++;
1841     }
1842   }
1843   delete [] yin0;
1844   delete [] yin1;
1845   return Float_t(naccept)/Float_t(nknots);
1846 }
1847
1848 Float_t  AliTPCcalibDButil::FilterTemperature(AliTPCSensorTempArray *tempArray, Double_t ymin, Double_t ymax, Double_t sigmaCut){
1849   //
1850   // Filter temperature array
1851   // tempArray    - array of temperatures         -
1852   // ymin         - minimal accepted temperature  - default 15
1853   // ymax         - maximal accepted temperature  - default 22
1854   // sigmaCut     - values filtered on interval median+-sigmaCut*rms - defaut 5
1855   // return value - fraction of filtered sensors
1856   const Double_t kMaxDy=0.1;
1857   Int_t nsensors=tempArray->NumSensors();
1858   if (nsensors==0) return 0.;
1859   Int_t naccept=0;
1860   for (Int_t isensor=0; isensor<nsensors; isensor++){
1861     AliDCSSensor *sensor = tempArray->GetSensorNum(isensor);
1862     if (!sensor) continue;
1863     //printf("%d\n",isensor);
1864     FilterSensor(sensor,ymin,ymax,kMaxDy, sigmaCut);
1865     if (sensor->GetFit()==0){
1866       //delete sensor;
1867       tempArray->RemoveSensorNum(isensor);
1868     }else{
1869       naccept++;
1870     }
1871   }
1872   return Float_t(naccept)/Float_t(nsensors);
1873 }
1874
1875
1876 void AliTPCcalibDButil::FilterCE(Double_t deltaT, Double_t cutAbs, Double_t cutSigma, TTreeSRedirector *pcstream){
1877   //
1878   // Filter CE data
1879   // Input parameters:
1880   //    deltaT   - smoothing window (in seconds)
1881   //    cutAbs   - max distance of the time info to the median (in time bins)
1882   //    cutSigma - max distance (in the RMS)
1883   //    pcstream - optional debug streamer to store original and filtered info
1884   // Hardwired parameters:
1885   //    kMinPoints =10;       // minimal number of points to define the CE
1886   //    kMinSectors=12;       // minimal number of sectors to define sideCE
1887   // Algorithm:
1888   // 0. Filter almost emty graphs (kMinPoints=10)
1889   // 1. calculate median and RMS per side
1890   // 2. Filter graphs - in respect with side medians
1891   //                  - cutAbs and cutDelta used
1892   // 3. Cut in respect wit the graph median - cutAbs and cutRMS used
1893   // 4. Calculate mean for A side and C side
1894   //
1895   const Int_t kMinPoints =10;       // minimal number of points to define the CE
1896   const Int_t kMinSectors=12;       // minimal number of sectors to define sideCE
1897   const Int_t kMinTime   =400;     // minimal arrival time of CE
1898   TObjArray *arrT=AliTPCcalibDB::Instance()->GetCErocTtime();
1899   Double_t medianY=0;
1900   TObjArray*  cearray =AliTPCcalibDB::Instance()->GetCEData();
1901   if (!cearray) return;
1902   Double_t tmin=-1;
1903   Double_t tmax=-1;
1904   //
1905   //
1906   AliTPCSensorTempArray *tempMapCE = (AliTPCSensorTempArray *)cearray->FindObject("TempMap");
1907   AliDCSSensor * cavernPressureCE  = (AliDCSSensor *) cearray->FindObject("CavernPressure");
1908   if ( tempMapCE && cavernPressureCE){
1909     //
1910     Bool_t isOK = FilterTemperature(tempMapCE)>0.1;
1911     FilterSensor(cavernPressureCE,960,1050,10, 5.);
1912     if (cavernPressureCE->GetFit()==0) isOK=kFALSE;
1913     if (isOK)  {
1914       // recalculate P/T correction map for time of the CE
1915       AliTPCCalibVdrift * driftCalib = new AliTPCCalibVdrift(tempMapCE,cavernPressureCE ,0);
1916       driftCalib->SetName("driftPTCE");
1917       driftCalib->SetTitle("driftPTCE");
1918       cearray->AddLast(driftCalib);
1919     }
1920   }
1921   //
1922   // 0. Filter almost emty graphs
1923   //
1924
1925   for (Int_t i=0; i<72;i++){
1926     TGraph *graph= (TGraph*)arrT->At(i);
1927     if (!graph) continue;
1928     if (graph->GetN()<kMinPoints){
1929       arrT->AddAt(0,i);
1930       delete graph;  // delete empty graph
1931       continue;
1932     }
1933     if (tmin<0) tmin = graph->GetX()[0];
1934     if (tmax<0) tmax = graph->GetX()[graph->GetN()-1];
1935     //
1936     if (tmin>graph->GetX()[0]) tmin=graph->GetX()[0];
1937     if (tmax<graph->GetX()[graph->GetN()-1]) tmax=graph->GetX()[graph->GetN()-1];
1938   }
1939   //
1940   // 1. calculate median and RMS per side
1941   //
1942   TArrayF arrA(100000), arrC(100000);
1943   Int_t nA=0, nC=0;
1944   Double_t medianA=0, medianC=0;
1945   Double_t rmsA=0, rmsC=0;
1946   for (Int_t isec=0; isec<72;isec++){
1947     TGraph *graph= (TGraph*)arrT->At(isec);
1948     if (!graph) continue;
1949     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1950       if (graph->GetY()[ipoint]<kMinTime) continue;
1951       if (nA>=arrA.fN) arrA.Set(nA*2);
1952       if (nC>=arrC.fN) arrC.Set(nC*2);
1953       if (isec%36<18)  arrA[nA++]= graph->GetY()[ipoint];
1954       if (isec%36>=18) arrC[nC++]= graph->GetY()[ipoint];
1955     }
1956   }
1957   if (nA>0){
1958     medianA=TMath::Median(nA,arrA.fArray);
1959     rmsA   =TMath::RMS(nA,arrA.fArray);
1960   }
1961   if (nC>0){
1962     medianC=TMath::Median(nC,arrC.fArray);
1963     rmsC   =TMath::RMS(nC,arrC.fArray);
1964   }
1965   //
1966   // 2. Filter graphs - in respect with side medians
1967   //
1968   TArrayD vecX(100000), vecY(100000);
1969   for (Int_t isec=0; isec<72;isec++){
1970     TGraph *graph= (TGraph*)arrT->At(isec);
1971     if (!graph) continue;
1972     Double_t median = (isec%36<18) ? medianA: medianC;
1973     Double_t rms    = (isec%36<18) ? rmsA:    rmsC;
1974     Int_t naccept=0;
1975     for (Int_t ipoint=kMinPoints-1; ipoint<graph->GetN();ipoint++){
1976       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutAbs) continue;
1977       if (TMath::Abs(graph->GetY()[ipoint]-median)>cutSigma*rms) continue;
1978       vecX[naccept]= graph->GetX()[ipoint];
1979       vecY[naccept]= graph->GetY()[ipoint];
1980       naccept++;
1981     }
1982     if (naccept<kMinPoints){
1983       arrT->AddAt(0,isec);
1984       delete graph;  // delete empty graph
1985       continue;
1986     }
1987     TGraph *graph2 = new TGraph(naccept, vecX.fArray, vecY.fArray);
1988     delete graph;
1989     arrT->AddAt(graph2,isec);
1990   }
1991   //
1992   // 3. Cut in respect wit the graph median
1993   //
1994   for (Int_t i=0; i<72;i++){
1995     TGraph *graph= (TGraph*)arrT->At(i);
1996     if (!graph) continue;
1997     //
1998     // filter in range
1999     //
2000     TGraph* graphTS0= FilterGraphMedianAbs(graph,cutAbs,medianY);
2001     if (!graphTS0) continue;
2002     if (graphTS0->GetN()<kMinPoints) {
2003       delete graphTS0;
2004       delete graph;
2005       arrT->AddAt(0,i);
2006       continue;
2007     }
2008     TGraph* graphTS= FilterGraphMedian(graphTS0,cutSigma,medianY);
2009     graphTS->Sort();
2010     AliTPCcalibDButil::SmoothGraph(graphTS,deltaT);
2011     if (pcstream){
2012       Int_t run = AliTPCcalibDB::Instance()->GetRun();
2013       (*pcstream)<<"filterCE"<<
2014         "run="<<run<<
2015         "isec="<<i<<
2016         "mY="<<medianY<<
2017         "graph.="<<graph<<
2018         "graphTS0.="<<graphTS0<<
2019         "graphTS.="<<graphTS<<
2020         "\n";
2021     }
2022     delete graphTS0;
2023     if (!graphTS) continue;
2024     arrT->AddAt(graphTS,i);
2025     delete graph;
2026   }
2027   //
2028   // Recalculate the mean time A side C side
2029   //
2030   TArrayF xA(200), yA(200), eA(200), xC(200),yC(200), eC(200);
2031   Int_t meanPoints=(nA+nC)/72;  // mean number of points
2032   for (Int_t itime=0; itime<200; itime++){
2033     nA=0, nC=0;
2034     Double_t time=tmin+(tmax-tmin)*Float_t(itime)/200.;
2035     for (Int_t i=0; i<72;i++){
2036       TGraph *graph= (TGraph*)arrT->At(i);
2037       if (!graph) continue;
2038       if (graph->GetN()<(meanPoints/4)) continue;
2039       if ( (i%36)<18 )  arrA[nA++]=graph->Eval(time);
2040       if ( (i%36)>=18 ) arrC[nC++]=graph->Eval(time);
2041     }
2042     xA[itime]=time;
2043     xC[itime]=time;
2044     yA[itime]=(nA>0)? TMath::Mean(nA,arrA.fArray):0;
2045     yC[itime]=(nC>0)? TMath::Mean(nC,arrC.fArray):0;
2046     eA[itime]=(nA>0)? TMath::RMS(nA,arrA.fArray):0;
2047     eC[itime]=(nC>0)? TMath::RMS(nC,arrC.fArray):0;
2048   }
2049   //
2050   Double_t rmsTA = TMath::RMS(200,yA.fArray)+TMath::Mean(200,eA.fArray);
2051   Double_t rmsTC = TMath::RMS(200,yC.fArray)+TMath::Mean(200,eC.fArray);
2052   if (pcstream){
2053     Int_t run = AliTPCcalibDB::Instance()->GetRun();
2054     (*pcstream)<<"filterAC"<<
2055       "run="<<run<<
2056       "nA="<<nA<<
2057       "nC="<<nC<<
2058       "rmsTA="<<rmsTA<<
2059       "rmsTC="<<rmsTC<<
2060       "\n";
2061   }
2062   //
2063   TGraphErrors *grA = new TGraphErrors(200,xA.fArray,yA.fArray,0, eA.fArray);
2064   TGraphErrors *grC = new TGraphErrors(200,xC.fArray,yC.fArray,0, eC.fArray);
2065   TGraph* graphTSA= FilterGraphMedian(grA,cutSigma,medianY);
2066   if (graphTSA&&graphTSA->GetN()) SmoothGraph(graphTSA,deltaT);
2067   TGraph* graphTSC= FilterGraphMedian(grC,cutSigma,medianY);
2068   if (graphTSC&&graphTSC->GetN()>0) SmoothGraph(graphTSC,deltaT);
2069   delete grA;
2070   delete grC;
2071   if (nA<kMinSectors) arrT->AddAt(0,72);
2072   else arrT->AddAt(graphTSA,72);
2073   if (nC<kMinSectors) arrT->AddAt(0,73);
2074   else arrT->AddAt(graphTSC,73);
2075 }
2076
2077
2078 void AliTPCcalibDButil::FilterTracks(Int_t run, Double_t cutSigma, TTreeSRedirector *pcstream){
2079   //
2080   // Filter Drift velocity measurement using the tracks
2081   // 0.  remove outlyers - error based
2082   //     cutSigma
2083   //
2084   //
2085   const Int_t kMinPoints=1;  // minimal number of points to define value
2086   TObjArray *arrT=AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2087   Double_t medianY=0;
2088   if (!arrT) return;
2089   for (Int_t i=0; i<arrT->GetEntries();i++){
2090     TGraphErrors *graph= (TGraphErrors*)arrT->At(i);
2091     if (!graph) continue;
2092     if (graph->GetN()<kMinPoints){
2093       delete graph;
2094       arrT->AddAt(0,i);
2095       continue;
2096     }
2097     TGraphErrors *graph2= FilterGraphMedianErr(graph,cutSigma,medianY);
2098     if (!graph2) {
2099       delete graph; arrT->AddAt(0,i); continue;
2100     }
2101     if (graph2->GetN()<1) {
2102       delete graph; arrT->AddAt(0,i); continue;
2103     }
2104     graph2->SetName(graph->GetName());
2105     graph2->SetTitle(graph->GetTitle());
2106     arrT->AddAt(graph2,i);
2107     if (pcstream){
2108       (*pcstream)<<"filterTracks"<<
2109         "run="<<run<<
2110         "isec="<<i<<
2111         "mY="<<medianY<<
2112         "graph.="<<graph<<
2113         "graph2.="<<graph2<<
2114         "\n";
2115     }
2116     delete graph;
2117   }
2118 }
2119
2120
2121
2122
2123
2124 Double_t AliTPCcalibDButil::GetLaserTime0(Int_t run, Int_t timeStamp, Int_t deltaT, Int_t side){
2125   //
2126   //
2127   // get laser time offset
2128   // median around timeStamp+-deltaT
2129   // QA - chi2 needed for later usage - to be added
2130   //    - currently cut on error
2131   //
2132   Int_t kMinPoints=1;
2133   Double_t kMinDelay=0.01;
2134   Double_t kMinDelayErr=0.0001;
2135
2136   TObjArray *array =AliTPCcalibDB::Instance()->GetTimeVdriftSplineRun(run);
2137   if (!array) return 0;
2138   TGraphErrors *tlaser=0;
2139   if (array){
2140     if (side==0) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_A");
2141     if (side==1) tlaser=(TGraphErrors*)array->FindObject("GRAPH_MEAN_DELAY_LASER_ALL_C");
2142   }
2143   if (!tlaser) return 0;
2144   Int_t npoints0= tlaser->GetN();
2145   if (npoints0==0) return 0;
2146   Double_t *xlaser = new Double_t[npoints0];
2147   Double_t *ylaser = new Double_t[npoints0];
2148   Int_t npoints=0;
2149   for (Int_t i=0;i<npoints0;i++){
2150     //printf("%d\n",i);
2151     if (tlaser->GetY()[i]<=kMinDelay) continue; // filter zeros
2152     if (tlaser->GetErrorY(i)>TMath::Abs(kMinDelayErr)) continue;
2153     xlaser[npoints]=tlaser->GetX()[npoints];
2154     ylaser[npoints]=tlaser->GetY()[npoints];
2155     npoints++;
2156   }
2157   //
2158   //
2159   Int_t index0=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp-deltaT))-1;
2160   Int_t index1=TMath::BinarySearch(npoints, xlaser, Double_t(timeStamp+deltaT))+1;
2161   //if (index1-index0 <kMinPoints) { index1+=kMinPoints; index0-=kMinPoints;}
2162   if (index0<0) index0=0;
2163   if (index1>=npoints-1) index1=npoints-1;
2164   if (index1-index0<kMinPoints) return 0;
2165   //
2166   //Double_t median = TMath::Median(index1-index0, &(ylaser[index0]));
2167     Double_t mean = TMath::Mean(index1-index0, &(ylaser[index0]));
2168   delete [] ylaser;
2169   delete [] xlaser;
2170   return mean;
2171 }
2172
2173
2174
2175
2176 void AliTPCcalibDButil::FilterGoofie(AliDCSSensorArray * goofieArray, Double_t deltaT, Double_t cutSigma, Double_t minVd, Double_t maxVd, TTreeSRedirector *pcstream){
2177   //
2178   // Filter Goofie data
2179   // goofieArray - points will be filtered
2180   // deltaT      - smmothing time window
2181   // cutSigma    - outler sigma cut in rms
2182   // minVn, maxVd- range absolute cut for variable vd/pt
2183   //             - to be tuned
2184   //
2185   // Ignore goofie if not enough points
2186   //
2187   const Int_t kMinPoints = 3;
2188   //
2189
2190   TGraph *graphvd = goofieArray->GetSensorNum(2)->GetGraph();
2191   TGraph *graphan = goofieArray->GetSensorNum(8)->GetGraph();
2192   TGraph *graphaf = goofieArray->GetSensorNum(9)->GetGraph();
2193   TGraph *graphpt = goofieArray->GetSensorNum(15)->GetGraph();
2194   if (!graphvd) return;
2195   if (graphvd->GetN()<kMinPoints){
2196     delete graphvd;
2197     goofieArray->GetSensorNum(2)->SetGraph(0);
2198     return;
2199   }
2200   //
2201   // 1. Caluclate medians of critical variables
2202   //    drift velcocity
2203   //    P/T
2204   //    area near peak
2205   //    area far  peak
2206   //
2207   Double_t medianpt=0;
2208   Double_t medianvd=0, sigmavd=0;
2209   Double_t medianan=0;
2210   Double_t medianaf=0;
2211   Int_t    entries=graphvd->GetN();
2212   Double_t yvdn[10000];
2213   Int_t nvd=0;
2214   //
2215   for (Int_t ipoint=0; ipoint<entries; ipoint++){
2216     if (graphpt->GetY()[ipoint]<=0.0000001) continue;
2217     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]<minVd) continue;
2218     if (graphvd->GetY()[ipoint]/graphpt->GetY()[ipoint]>maxVd) continue;
2219     yvdn[nvd++]=graphvd->GetY()[ipoint];
2220   }
2221   if (nvd<kMinPoints){
2222     delete graphvd;
2223     goofieArray->GetSensorNum(2)->SetGraph(0);
2224     return;
2225   }
2226   //
2227   Int_t nuni = TMath::Min(TMath::Nint(nvd*0.4+2), nvd-1);
2228   if (nuni>=kMinPoints){
2229     AliMathBase::EvaluateUni(nvd, yvdn, medianvd,sigmavd,nuni);
2230   }else{
2231     medianvd = TMath::Median(nvd, yvdn);
2232   }
2233
2234   TGraph * graphpt0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphpt,10,medianpt);
2235   TGraph * graphpt1 = AliTPCcalibDButil::FilterGraphMedian(graphpt0,2,medianpt);
2236   TGraph * graphan0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphan,10,medianan);
2237   TGraph * graphan1 = AliTPCcalibDButil::FilterGraphMedian(graphan0,2,medianan);
2238   TGraph * graphaf0 = AliTPCcalibDButil::FilterGraphMedianAbs(graphaf,10,medianaf);
2239   TGraph * graphaf1 = AliTPCcalibDButil::FilterGraphMedian(graphaf0,2,medianaf);
2240   delete graphpt0;
2241   delete graphpt1;
2242   delete graphan0;
2243   delete graphan1;
2244   delete graphaf0;
2245   delete graphaf1;
2246   //
2247   // 2. Make outlyer graph
2248   //
2249   Int_t nOK=0;
2250   TGraph graphOut(*graphvd);
2251   for (Int_t i=0; i<entries;i++){
2252     //
2253     Bool_t isOut=kFALSE;
2254     if (graphpt->GetY()[i]<=0.0000001) {  graphOut.GetY()[i]=1; continue;}
2255     if (graphvd->GetY()[i]/graphpt->GetY()[i]<minVd || graphvd->GetY()[i]/graphpt->GetY()[i]>maxVd) {  graphOut.GetY()[i]=1; continue;}
2256
2257     if (TMath::Abs((graphvd->GetY()[i]/graphpt->GetY()[i])/medianvd-1.)<0.05)
2258       isOut|=kTRUE;
2259     if (TMath::Abs(graphpt->GetY()[i]/medianpt-1.)>0.02) isOut|=kTRUE;
2260     if (TMath::Abs(graphan->GetY()[i]/medianan-1.)>0.2) isOut|=kTRUE;
2261     if (TMath::Abs(graphaf->GetY()[i]/medianaf-1.)>0.2) isOut|=kTRUE;
2262     graphOut.GetY()[i]= (isOut)?1:0;
2263     if (!isOut) nOK++;
2264   }
2265   if (nOK<kMinPoints) {
2266     delete graphvd;
2267     goofieArray->GetSensorNum(2)->SetGraph(0);
2268     return;
2269   }
2270   //
2271   // 3. Filter out outlyers - and smooth
2272   //
2273   TVectorF vmedianArray(goofieArray->NumSensors());
2274   TVectorF vrmsArray(goofieArray->NumSensors());
2275   Double_t xnew[10000];
2276   Double_t ynew[10000];
2277   TObjArray junk;
2278   junk.SetOwner(kTRUE);
2279   Bool_t isOK=kTRUE;
2280   //
2281   //
2282   for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
2283     isOK=kTRUE;
2284     AliDCSSensor *sensor = goofieArray->GetSensorNum(isensor);
2285     TGraph *graphOld=0, *graphNew=0, * graphNew0=0,*graphNew1=0,*graphNew2=0;
2286     //
2287     if (!sensor) continue;
2288     graphOld = sensor->GetGraph();
2289     if (graphOld) {
2290       sensor->SetGraph(0);
2291       Int_t nused=0;
2292       for (Int_t i=0;i<entries;i++){
2293         if (graphOut.GetY()[i]>0.5) continue;
2294         xnew[nused]=graphOld->GetX()[i];
2295         ynew[nused]=graphOld->GetY()[i];
2296         nused++;
2297       }
2298       graphNew = new TGraph(nused,xnew,ynew);
2299       junk.AddLast(graphNew);
2300       junk.AddLast(graphOld);
2301       Double_t median=0;
2302       graphNew0  = AliTPCcalibDButil::FilterGraphMedian(graphNew,cutSigma,median);
2303       if (graphNew0!=0){
2304         junk.AddLast(graphNew0);
2305         graphNew1  = AliTPCcalibDButil::FilterGraphMedian(graphNew0,cutSigma,median);
2306         if (graphNew1!=0){
2307           junk.AddLast(graphNew1);
2308           graphNew2  = AliTPCcalibDButil::FilterGraphMedian(graphNew1,cutSigma,median);
2309           if (graphNew2!=0) {
2310             vrmsArray[isensor]   =TMath::RMS(graphNew2->GetN(),graphNew2->GetY());
2311             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2312             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2313             AliTPCcalibDButil::SmoothGraph(graphNew2,deltaT);
2314             printf("%d\t%f\t%f\n",isensor, median,vrmsArray[isensor]);
2315             vmedianArray[isensor]=median;
2316             //
2317           }
2318         }
2319       }
2320     }
2321     if (!graphOld) {  isOK=kFALSE; graphOld =&graphOut;}
2322     if (!graphNew0) { isOK=kFALSE; graphNew0=graphOld;}
2323     if (!graphNew1) { isOK=kFALSE; graphNew1=graphOld;}
2324     if (!graphNew2) { isOK=kFALSE; graphNew2=graphOld;}
2325     (*pcstream)<<"goofieA"<<
2326       Form("isOK_%d.=",isensor)<<isOK<<
2327       Form("s_%d.=",isensor)<<sensor<<
2328       Form("gr_%d.=",isensor)<<graphOld<<
2329       Form("gr0_%d.=",isensor)<<graphNew0<<
2330       Form("gr1_%d.=",isensor)<<graphNew1<<
2331       Form("gr2_%d.=",isensor)<<graphNew2;
2332     if (isOK) sensor->SetGraph(graphNew2);
2333   }
2334   (*pcstream)<<"goofieA"<<
2335     "vmed.="<<&vmedianArray<<
2336     "vrms.="<<&vrmsArray<<
2337     "\n";
2338   junk.Delete();   // delete temoprary graphs
2339
2340 }
2341
2342
2343
2344
2345
2346