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
  34 #include <AliDCSSensorArray.h>
  35 #include <AliDCSSensor.h>
  36 #include "AliTPCcalibDB.h"
  37 #include "AliTPCCalPad.h"
  38 #include "AliTPCCalROC.h"
  39 #include "AliTPCROC.h"
  40 #include "AliTPCmapper.h"
  41 #include "AliTPCParam.h"
  42 #include "AliTPCCalibRaw.h"
  43
  44 #include "AliTPCcalibDButil.h"
  45
  46 ClassImp(AliTPCcalibDButil)
  47 AliTPCcalibDButil::AliTPCcalibDButil() :
  48   TObject(),
  49   fCalibDB(AliTPCcalibDB::Instance()),
  50   fPadNoise(0x0),
  51   fPedestals(0x0),
  52   fPulserTmean(0x0),
  53   fPulserTrms(0x0),
  54   fPulserQmean(0x0),
  55   fPulserOutlier(new AliTPCCalPad("PulserOutliers","PulserOutliers")),
  56   fCETmean(0x0),
  57   fCETrms(0x0),
  58   fCEQmean(0x0),
  59   fALTROMasked(0x0),
  60   fCalibRaw(0x0),
  61   fRefPadNoise(0x0),
  62   fRefPedestals(0x0),
  63   fRefPulserTmean(0x0),
  64   fRefPulserTrms(0x0),
  65   fRefPulserQmean(0x0),
  66   fRefPulserOutlier(new AliTPCCalPad("RefPulserOutliers","RefPulserOutliers")),
  67   fRefCETmean(0x0),
  68   fRefCETrms(0x0),
  69   fRefCEQmean(0x0),
  70   fRefALTROMasked(0x0),
  71   fRefCalibRaw(0x0),
  72   fGoofieArray(0x0),
  73   fMapper(new AliTPCmapper(0x0)),
  74   fNpulserOutliers(-1),
  75   fIrocTimeOffset(0),
  76   fCETmaxLimitAbs(1.5),
  77   fPulTmaxLimitAbs(1.5),
  78   fPulQmaxLimitAbs(5),
  79   fPulQminLimit(11)
  80 {
  81   //
  82   // Default ctor
  83   //
  84 }
  85 //_____________________________________________________________________________________
  86 AliTPCcalibDButil::~AliTPCcalibDButil()
  87 {
  88   //
  89   // dtor
  90   //
  91   delete fPulserOutlier;
  92   delete fRefPulserOutlier;
  93   delete fMapper;
  94   if (fRefPadNoise) delete fRefPadNoise;
  95   if (fRefPedestals) delete fRefPedestals;
  96   if (fRefPulserTmean) delete fRefPulserTmean;
  97   if (fRefPulserTrms) delete fRefPulserTrms;
  98   if (fRefPulserQmean) delete fRefPulserQmean;
  99   if (fRefCETmean) delete fRefCETmean;
 100   if (fRefCETrms) delete fRefCETrms;
 101   if (fRefCEQmean) delete fRefCEQmean;
 102   if (fRefALTROMasked) delete fRefALTROMasked;
 103   if (fRefCalibRaw) delete fRefCalibRaw;
 104
 105 }
 106 //_____________________________________________________________________________________
 107 void AliTPCcalibDButil::UpdateFromCalibDB()
 108 {
 109   //
 110   // Update pointers from calibDB
 111   //
 112   fPadNoise=fCalibDB->GetPadNoise();
 113   fPedestals=fCalibDB->GetPedestals();
 114   fPulserTmean=fCalibDB->GetPulserTmean();
 115   fPulserTrms=fCalibDB->GetPulserTrms();
 116   fPulserQmean=fCalibDB->GetPulserQmean();
 117   fCETmean=fCalibDB->GetCETmean();
 118   fCETrms=fCalibDB->GetCETrms();
 119   fCEQmean=fCalibDB->GetCEQmean();
 120   fALTROMasked=fCalibDB->GetALTROMasked();
 121   fGoofieArray=fCalibDB->GetGoofieSensors(fCalibDB->GetRun());
 122   fCalibRaw=fCalibDB->GetCalibRaw();
 123   UpdatePulserOutlierMap();
 124 }
 125 //_____________________________________________________________________________________
 126 void AliTPCcalibDButil::ProcessCEdata(const char* fitFormula, TVectorD &fitResultsA, TVectorD &fitResultsC,
 127                                       Int_t &noutliersCE, AliTPCCalPad *outCE)
 128 {
 129   //
 130   // Process the CE data for this run
 131   // the return TVectorD arrays contian the results of the fit
 132   // noutliersCE contains the number of pads marked as outliers,
 133   //   not including masked and edge pads
 134   //
 135
 136   //retrieve CE and ALTRO data
 137   if (!fCETmean){
 138     TString fitString(fitFormula);
 139     fitString.ReplaceAll("++","#");
 140     Int_t ndim=fitString.CountChar('#')+2;
 141     fitResultsA.ResizeTo(ndim);
 142     fitResultsC.ResizeTo(ndim);
 143     fitResultsA.Zero();
 144     fitResultsC.Zero();
 145     noutliersCE=-1;
 146     return;
 147   }
 148   noutliersCE=0;
 149   //create outlier map
 150   AliTPCCalPad *out=0;
 151   if (outCE) out=outCE;
 152   else out=new AliTPCCalPad("outCE","outCE");
 153   AliTPCCalROC *rocMasked=0x0;
 154   //loop over all channels
 155   for (UInt_t iroc=0;iroc<fCETmean->kNsec;++iroc){
 156     AliTPCCalROC *rocData=fCETmean->GetCalROC(iroc);
 157     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(iroc);
 158     AliTPCCalROC *rocOut=out->GetCalROC(iroc);
 159     if (!rocData) {
 160       noutliersCE+=AliTPCROC::Instance()->GetNChannels(iroc);
 161       rocOut->Add(1.);
 162       continue;
 163     }
 164     //add time offset to IROCs
 165     if (iroc<AliTPCROC::Instance()->GetNInnerSector())
 166       rocData->Add(fIrocTimeOffset);
 167     //select outliers
 168     UInt_t nrows=rocData->GetNrows();
 169     for (UInt_t irow=0;irow<nrows;++irow){
 170       UInt_t npads=rocData->GetNPads(irow);
 171       for (UInt_t ipad=0;ipad<npads;++ipad){
 172         rocOut->SetValue(irow,ipad,0);
 173         //exclude masked pads
 174         if (rocMasked && rocMasked->GetValue(irow,ipad)) {
 175           rocOut->SetValue(irow,ipad,1);
 176           continue;
 177         }
 178         //exclude edge pads
 179         if (ipad==0||ipad==npads-1) rocOut->SetValue(irow,ipad,1);
 180         Float_t valTmean=rocData->GetValue(irow,ipad);
 181         //exclude values that are exactly 0
 182         if (valTmean==0) {
 183           rocOut->SetValue(irow,ipad,1);
 184           ++noutliersCE;
 185         }
 186         // exclude channels with too large variations
 187         if (TMath::Abs(valTmean)>fCETmaxLimitAbs) {
 188           rocOut->SetValue(irow,ipad,1);
 189           ++noutliersCE;
 190         }
 191       }
 192     }
 193   }
 194   //perform fit
 195   TMatrixD dummy;
 196   Float_t chi2A,chi2C;
 197   fCETmean->GlobalSidesFit(out,fitFormula,fitResultsA,fitResultsC,dummy,dummy,chi2A,chi2C);
 198   if (!outCE) delete out;
 199 }
 200 //_____________________________________________________________________________________
 201 void AliTPCcalibDButil::ProcessCEgraphs(TVectorD &vecTEntries, TVectorD &vecTMean, TVectorD &vecTRMS, TVectorD &vecTMedian,
 202                      TVectorD &vecQEntries, TVectorD &vecQMean, TVectorD &vecQRMS, TVectorD &vecQMedian,
 203                      Float_t &driftTimeA, Float_t &driftTimeC )
 204 {
 205   //
 206   // Calculate statistical information from the CE graphs for drift time and charge
 207   //
 208
 209   //reset arrays
 210   vecTEntries.ResizeTo(72);
 211   vecTMean.ResizeTo(72);
 212   vecTRMS.ResizeTo(72);
 213   vecTMedian.ResizeTo(72);
 214   vecQEntries.ResizeTo(72);
 215   vecQMean.ResizeTo(72);
 216   vecQRMS.ResizeTo(72);
 217   vecQMedian.ResizeTo(72);
 218   vecTEntries.Zero();
 219   vecTMean.Zero();
 220   vecTRMS.Zero();
 221   vecTMedian.Zero();
 222   vecQEntries.Zero();
 223   vecQMean.Zero();
 224   vecQRMS.Zero();
 225   vecQMedian.Zero();
 226   driftTimeA=0;
 227   driftTimeC=0;
 228   TObjArray *arrT=fCalibDB->GetCErocTtime();
 229   TObjArray *arrQ=fCalibDB->GetCErocQtime();
 230   if (arrT){
 231     for (Int_t isec=0;isec<74;++isec){
 232       TGraph *gr=(TGraph*)arrT->At(isec);
 233       if (!gr) continue;
 234       TVectorD values;
 235       Int_t npoints = gr->GetN();
 236       values.ResizeTo(npoints);
 237       Int_t nused =0;
 238       //skip first points, theres always some problems with finding the CE position
 239       for (Int_t ipoint=4; ipoint<npoints; ipoint++){
 240         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1020 ){
 241           values[nused]=gr->GetY()[ipoint];
 242           nused++;
 243         }
 244       }
 245       //
 246       if (isec<72) vecTEntries[isec]= nused;
 247       if (nused>1){
 248         if (isec<72){
 249           vecTMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 250           vecTMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 251           vecTRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 252         } else if (isec==72){
 253           driftTimeA=TMath::Median(nused,values.GetMatrixArray());
 254         } else if (isec==73){
 255           driftTimeC=TMath::Median(nused,values.GetMatrixArray());
 256         }
 257       }
 258     }
 259   }
 260   if (arrQ){
 261     for (Int_t isec=0;isec<arrQ->GetEntriesFast();++isec){
 262       TGraph *gr=(TGraph*)arrQ->At(isec);
 263       if (!gr) continue;
 264       TVectorD values;
 265       Int_t npoints = gr->GetN();
 266       values.ResizeTo(npoints);
 267       Int_t nused =0;
 268       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 269         if (gr->GetY()[ipoint]>500 && gr->GetY()[ipoint]<1000 ){
 270           values[nused]=gr->GetY()[ipoint];
 271           nused++;
 272         }
 273       }
 274       //
 275       vecQEntries[isec]= nused;
 276       if (nused>1){
 277         vecQMedian[isec] = TMath::Median(nused,values.GetMatrixArray());
 278         vecQMean[isec]   = TMath::Mean(nused,values.GetMatrixArray());
 279         vecQRMS[isec]    = TMath::RMS(nused,values.GetMatrixArray());
 280       }
 281     }
 282   }
 283 }
 284
 285 //_____________________________________________________________________________________
 286 void AliTPCcalibDButil::ProcessNoiseData(TVectorD &vNoiseMean, TVectorD &vNoiseMeanSenRegions,
 287                       TVectorD &vNoiseRMS, TVectorD &vNoiseRMSSenRegions,
 288                       Int_t &nonMaskedZero)
 289 {
 290   //
 291   // process noise data
 292   // vNoiseMean/RMS contains the Mean/RMS noise of the complete TPC [0], IROCs only [1],
 293   //    OROCs small pads [2] and OROCs large pads [3]
 294   // vNoiseMean/RMSsenRegions constains the same information, but only for the sensitive regions (edge pads, corners, IROC spot)
 295   // nonMaskedZero contains the number of pads which show zero noise and were not masked. This might indicate an error
 296   //
 297
 298   //set proper size and reset
 299   const UInt_t infoSize=4;
 300   vNoiseMean.ResizeTo(infoSize);
 301   vNoiseMeanSenRegions.ResizeTo(infoSize);
 302   vNoiseRMS.ResizeTo(infoSize);
 303   vNoiseRMSSenRegions.ResizeTo(infoSize);
 304   vNoiseMean.Zero();
 305   vNoiseMeanSenRegions.Zero();
 306   vNoiseRMS.Zero();
 307   vNoiseRMSSenRegions.Zero();
 308   nonMaskedZero=0;
 309   //counters
 310   TVectorD c(infoSize);
 311   TVectorD cs(infoSize);
 312   //tpc parameters
 313   AliTPCParam par;
 314   par.Update();
 315   //retrieve noise and ALTRO data
 316   if (!fPadNoise) return;
 317   AliTPCCalROC *rocMasked=0x0;
 318   //create IROC, OROC1, OROC2 and sensitive region masks
 319   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 320     AliTPCCalROC *noiseROC=fPadNoise->GetCalROC(isec);
 321     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 322     UInt_t nrows=noiseROC->GetNrows();
 323     for (UInt_t irow=0;irow<nrows;++irow){
 324       UInt_t npads=noiseROC->GetNPads(irow);
 325       for (UInt_t ipad=0;ipad<npads;++ipad){
 326         //don't use masked channels;
 327         if (rocMasked && rocMasked->GetValue(irow,ipad)) continue;
 328         Float_t noiseVal=noiseROC->GetValue(irow,ipad);
 329         //check if noise==0
 330         if (noiseVal==0) {
 331           ++nonMaskedZero;
 332           continue;
 333         }
 334         //check for nan
 335         if ( !(noiseVal<10000000) ){
 336           printf ("Warning: nan detected in (sec,row,pad - val): %02d,%02d,%03d - %.1f\n",isec,irow,ipad,noiseVal);
 337           continue;
 338         }
 339         Int_t cpad=(Int_t)ipad-(Int_t)npads/2;
 340         Int_t masksen=1; // sensitive pards are not masked (0)
 341         if (ipad<2||npads-ipad-1<2) masksen=0; //don't mask edge pads (sensitive)
 342         if (isec<AliTPCROC::Instance()->GetNInnerSector()){
 343           //IROCs
 344           if (irow>19&&irow<46){
 345             if (TMath::Abs(cpad)<7) masksen=0; //IROC spot
 346           }
 347           Int_t type=1;
 348           vNoiseMean[type]+=noiseVal;
 349           vNoiseRMS[type]+=noiseVal*noiseVal;
 350           ++c[type];
 351           if (!masksen){
 352             vNoiseMeanSenRegions[type]+=noiseVal;
 353             vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 354             ++cs[type];
 355           }
 356         } else {
 357           //OROCs
 358           //define sensive regions
 359           if ((nrows-irow-1)<3) masksen=0; //last three rows in OROCs are sensitive
 360           if ( irow>75 ){
 361             Int_t padEdge=(Int_t)TMath::Min(ipad,npads-ipad);
 362             if (padEdge<((((Int_t)irow-76)/4+1))*2) masksen=0; //OROC outer corners are sensitive
 363           }
 364           if ((Int_t)irow<par.GetNRowUp1()){
 365             //OROC1
 366             Int_t type=2;
 367             vNoiseMean[type]+=noiseVal;
 368             vNoiseRMS[type]+=noiseVal*noiseVal;
 369             ++c[type];
 370             if (!masksen){
 371               vNoiseMeanSenRegions[type]+=noiseVal;
 372               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 373               ++cs[type];
 374             }
 375           }else{
 376             //OROC2
 377             Int_t type=3;
 378             vNoiseMean[type]+=noiseVal;
 379             vNoiseRMS[type]+=noiseVal*noiseVal;
 380             ++c[type];
 381             if (!masksen){
 382               vNoiseMeanSenRegions[type]+=noiseVal;
 383               vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 384               ++cs[type];
 385             }
 386           }
 387         }
 388         //whole tpc
 389         Int_t type=0;
 390         vNoiseMean[type]+=noiseVal;
 391         vNoiseRMS[type]+=noiseVal*noiseVal;
 392         ++c[type];
 393         if (!masksen){
 394           vNoiseMeanSenRegions[type]+=noiseVal;
 395           vNoiseRMSSenRegions[type]+=noiseVal*noiseVal;
 396           ++cs[type];
 397         }
 398       }//end loop pads
 399     }//end loop rows
 400   }//end loop sectors (rocs)
 401
 402   //calculate mean and RMS
 403   const Double_t verySmall=0.0000000001;
 404   for (UInt_t i=0;i<infoSize;++i){
 405     Double_t mean=0;
 406     Double_t rms=0;
 407     Double_t meanSen=0;
 408     Double_t rmsSen=0;
 409
 410     if (c[i]>verySmall){
 411 //       printf ("i: %d - m: %.3f, c: %.0f, r: %.3f\n",i,vNoiseMean[i],c[i],vNoiseRMS[i]);
 412       mean=vNoiseMean[i]/c[i];
 413       rms=vNoiseRMS[i];
 414       rms=TMath::Sqrt(TMath::Abs(rms/c[i]-mean*mean));
 415     }
 416     vNoiseMean[i]=mean;
 417     vNoiseRMS[i]=rms;
 418
 419     if (cs[i]>verySmall){
 420       meanSen=vNoiseMeanSenRegions[i]/cs[i];
 421       rmsSen=vNoiseRMSSenRegions[i];
 422       rmsSen=TMath::Sqrt(TMath::Abs(rmsSen/cs[i]-meanSen*meanSen));
 423     }
 424     vNoiseMeanSenRegions[i]=meanSen;
 425     vNoiseRMSSenRegions[i]=rmsSen;
 426   }
 427 }
 428
 429 //_____________________________________________________________________________________
 430 void AliTPCcalibDButil::ProcessPulser(TVectorD &vMeanTime)
 431 {
 432   //
 433   // Process the Pulser information
 434   // vMeanTime:     pulser mean time position in IROC-A, IROC-C, OROC-A, OROC-C
 435   //
 436
 437   const UInt_t infoSize=4;
 438   //reset counters to error number
 439   vMeanTime.ResizeTo(infoSize);
 440   vMeanTime.Zero();
 441   //counter
 442   TVectorD c(infoSize);
 443   //retrieve pulser and ALTRO data
 444   if (!fPulserTmean) return;
 445   //
 446   //get Outliers
 447   AliTPCCalROC *rocOut=0x0;
 448   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 449     AliTPCCalROC *tmeanROC=fPulserTmean->GetCalROC(isec);
 450     if (!tmeanROC) continue;
 451     rocOut=fPulserOutlier->GetCalROC(isec);
 452     UInt_t nchannels=tmeanROC->GetNchannels();
 453     for (UInt_t ichannel=0;ichannel<nchannels;++ichannel){
 454       if (rocOut && rocOut->GetValue(ichannel)) continue;
 455       Float_t val=tmeanROC->GetValue(ichannel);
 456       Int_t type=isec/18;
 457       vMeanTime[type]+=val;
 458       ++c[type];
 459     }
 460   }
 461   //calculate mean
 462   for (UInt_t itype=0; itype<infoSize; ++itype){
 463     if (c[itype]>0) vMeanTime[itype]/=c[itype];
 464     else vMeanTime[itype]=0;
 465   }
 466 }
 467 //_____________________________________________________________________________________
 468 void AliTPCcalibDButil::ProcessALTROConfig(Int_t &nMasked)
 469 {
 470   //
 471   // Get Values from ALTRO configuration data
 472   //
 473   nMasked=-1;
 474   if (!fALTROMasked) return;
 475   nMasked=0;
 476   for (Int_t isec=0;isec<fALTROMasked->kNsec; ++isec){
 477     AliTPCCalROC *rocMasked=fALTROMasked->GetCalROC(isec);
 478     for (UInt_t ichannel=0; ichannel<rocMasked->GetNchannels();++ichannel){
 479       if (rocMasked->GetValue(ichannel)) ++nMasked;
 480     }
 481   }
 482 }
 483 //_____________________________________________________________________________________
 484 void AliTPCcalibDButil::ProcessGoofie(TVectorD & vecEntries, TVectorD & vecMedian, TVectorD &vecMean, TVectorD &vecRMS)
 485 {
 486   //
 487   // Proces Goofie values, return statistical information of the currently set goofieArray
 488   // The meaning of the entries are given below
 489   /*
 490   1       TPC_ANODE_I_A00_STAT
 491   2       TPC_DVM_CO2
 492   3       TPC_DVM_DriftVelocity
 493   4       TPC_DVM_FCageHV
 494   5       TPC_DVM_GainFar
 495   6       TPC_DVM_GainNear
 496   7       TPC_DVM_N2
 497   8       TPC_DVM_NumberOfSparks
 498   9       TPC_DVM_PeakAreaFar
 499   10      TPC_DVM_PeakAreaNear
 500   11      TPC_DVM_PeakPosFar
 501   12      TPC_DVM_PeakPosNear
 502   13      TPC_DVM_PickupHV
 503   14      TPC_DVM_Pressure
 504   15      TPC_DVM_T1_Over_P
 505   16      TPC_DVM_T2_Over_P
 506   17      TPC_DVM_T_Over_P
 507   18      TPC_DVM_TemperatureS1
 508    */
 509   if (!fGoofieArray){
 510     Int_t nsensors=19;
 511     vecEntries.ResizeTo(nsensors);
 512     vecMedian.ResizeTo(nsensors);
 513     vecMean.ResizeTo(nsensors);
 514     vecRMS.ResizeTo(nsensors);
 515     vecEntries.Zero();
 516     vecMedian.Zero();
 517     vecMean.Zero();
 518     vecRMS.Zero();
 519     return;
 520   }
 521   Double_t kEpsilon=0.0000000001;
 522   Double_t kBig=100000000000.;
 523   Int_t nsensors = fGoofieArray->NumSensors();
 524   vecEntries.ResizeTo(nsensors);
 525   vecMedian.ResizeTo(nsensors);
 526   vecMean.ResizeTo(nsensors);
 527   vecRMS.ResizeTo(nsensors);
 528   TVectorF values;
 529   for (Int_t isensor=0; isensor<fGoofieArray->NumSensors();isensor++){
 530     AliDCSSensor *gsensor = fGoofieArray->GetSensor(isensor);
 531     if (gsensor &&  gsensor->GetGraph()){
 532       Int_t npoints = gsensor->GetGraph()->GetN();
 533       // filter zeroes
 534       values.ResizeTo(npoints);
 535       Int_t nused =0;
 536       for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 537         if (TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])>kEpsilon &&
 538             TMath::Abs(gsensor->GetGraph()->GetY()[ipoint])<kBig ){
 539               values[nused]=gsensor->GetGraph()->GetY()[ipoint];
 540               nused++;
 541             }
 542       }
 543       //
 544       vecEntries[isensor]= nused;
 545       if (nused>1){
 546         vecMedian[isensor] = TMath::Median(nused,values.GetMatrixArray());
 547         vecMean[isensor]   = TMath::Mean(nused,values.GetMatrixArray());
 548         vecRMS[isensor]    = TMath::RMS(nused,values.GetMatrixArray());
 549       }
 550     }
 551   }
 552 }
 553 //_____________________________________________________________________________________
 554 void AliTPCcalibDButil::ProcessPedestalVariations(TVectorF &pedestalDeviations)
 555 {
 556   //
 557   // check the variations of the pedestal data to the reference pedestal data
 558   // thresholds are 0.5, 1.0, 1.5 and 2 timebins respectively.
 559   //
 560   const Int_t npar=4;
 561   TVectorF vThres(npar); //thresholds
 562   Int_t nActive=0;       //number of active channels
 563
 564   //reset and set thresholds
 565   pedestalDeviations.ResizeTo(npar);
 566   for (Int_t i=0;i<npar;++i){
 567     pedestalDeviations.GetMatrixArray()[i]=0;
 568     vThres.GetMatrixArray()[i]=(i+1)*.5;
 569   }
 570   //check all needed data is available
 571   if (!fRefPedestals || !fPedestals || !fALTROMasked || !fRefALTROMasked) return;
 572   //loop over all channels
 573   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 574     AliTPCCalROC *pROC=fPedestals->GetCalROC(isec);
 575     AliTPCCalROC *pRefROC=fRefPedestals->GetCalROC(isec);
 576     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 577     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 578     UInt_t nrows=mROC->GetNrows();
 579     for (UInt_t irow=0;irow<nrows;++irow){
 580       UInt_t npads=mROC->GetNPads(irow);
 581       for (UInt_t ipad=0;ipad<npads;++ipad){
 582         //don't use masked channels;
 583         if (mROC   ->GetValue(irow,ipad)) continue;
 584         if (mRefROC->GetValue(irow,ipad)) continue;
 585         Float_t deviation=TMath::Abs(pROC->GetValue(irow,ipad)-pRefROC->GetValue(irow,ipad));
 586         for (Int_t i=0;i<npar;++i){
 587           if (deviation>vThres[i])
 588             ++pedestalDeviations.GetMatrixArray()[i];
 589         }
 590         ++nActive;
 591       }//end ipad
 592     }//ind irow
 593   }//end isec
 594   if (nActive>0){
 595     for (Int_t i=0;i<npar;++i){
 596       pedestalDeviations.GetMatrixArray()[i]/=nActive;
 597     }
 598   }
 599 }
 600 //_____________________________________________________________________________________
 601 void AliTPCcalibDButil::ProcessNoiseVariations(TVectorF &noiseDeviations)
 602 {
 603   //
 604   // check the variations of the noise data to the reference noise data
 605   // thresholds are 5, 10, 15 and 20 percent respectively.
 606   //
 607   const Int_t npar=4;
 608   TVectorF vThres(npar); //thresholds
 609   Int_t nActive=0;       //number of active channels
 610
 611   //reset and set thresholds
 612   noiseDeviations.ResizeTo(npar);
 613   for (Int_t i=0;i<npar;++i){
 614     noiseDeviations.GetMatrixArray()[i]=0;
 615     vThres.GetMatrixArray()[i]=(i+1)*.05;
 616   }
 617   //check all needed data is available
 618   if (!fRefPadNoise || !fPadNoise || !fALTROMasked || !fRefALTROMasked) return;
 619   //loop over all channels
 620   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 621     AliTPCCalROC *nROC=fPadNoise->GetCalROC(isec);
 622     AliTPCCalROC *nRefROC=fRefPadNoise->GetCalROC(isec);
 623     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 624     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 625     UInt_t nrows=mROC->GetNrows();
 626     for (UInt_t irow=0;irow<nrows;++irow){
 627       UInt_t npads=mROC->GetNPads(irow);
 628       for (UInt_t ipad=0;ipad<npads;++ipad){
 629         //don't use masked channels;
 630         if (mROC   ->GetValue(irow,ipad)) continue;
 631         if (mRefROC->GetValue(irow,ipad)) continue;
 632         Float_t deviation=(nROC->GetValue(irow,ipad)/nRefROC->GetValue(irow,ipad))-1;
 633         for (Int_t i=0;i<npar;++i){
 634           if (deviation>vThres[i])
 635             ++noiseDeviations.GetMatrixArray()[i];
 636         }
 637         ++nActive;
 638       }//end ipad
 639     }//ind irow
 640   }//end isec
 641   if (nActive>0){
 642     for (Int_t i=0;i<npar;++i){
 643       noiseDeviations.GetMatrixArray()[i]/=nActive;
 644     }
 645   }
 646 }
 647 //_____________________________________________________________________________________
 648 void AliTPCcalibDButil::ProcessPulserVariations(TVectorF &pulserQdeviations, Float_t &varQMean,
 649                                                 Int_t &npadsOutOneTB, Int_t &npadsOffAdd)
 650 {
 651   //
 652   // check the variations of the pulserQmean data to the reference pulserQmean data: pulserQdeviations
 653   // thresholds are .5, 1, 5 and 10 percent respectively.
 654   //
 655   //
 656   const Int_t npar=4;
 657   TVectorF vThres(npar); //thresholds
 658   Int_t nActive=0;       //number of active channels
 659
 660   //reset and set thresholds
 661   pulserQdeviations.ResizeTo(npar);
 662   for (Int_t i=0;i<npar;++i){
 663     pulserQdeviations.GetMatrixArray()[i]=0;
 664   }
 665   npadsOutOneTB=0;
 666   npadsOffAdd=0;
 667   varQMean=0;
 668   vThres.GetMatrixArray()[0]=.005;
 669   vThres.GetMatrixArray()[1]=.01;
 670   vThres.GetMatrixArray()[2]=.05;
 671   vThres.GetMatrixArray()[3]=.1;
 672   //check all needed data is available
 673   if (!fRefPulserTmean || !fPulserTmean || !fPulserQmean || !fRefPulserQmean || !fALTROMasked || !fRefALTROMasked) return;
 674   //
 675   UpdateRefPulserOutlierMap();
 676   //loop over all channels
 677   for (UInt_t isec=0;isec<(UInt_t)AliTPCCalPad::kNsec;++isec){
 678     AliTPCCalROC *pqROC=fPulserQmean->GetCalROC(isec);
 679     AliTPCCalROC *pqRefROC=fRefPulserQmean->GetCalROC(isec);
 680     AliTPCCalROC *ptROC=fPulserTmean->GetCalROC(isec);
 681     //AliTPCCalROC *ptRefROC=fRefPulserTmean->GetCalROC(isec);
 682     AliTPCCalROC *mROC=fALTROMasked->GetCalROC(isec);
 683     AliTPCCalROC *mRefROC=fRefALTROMasked->GetCalROC(isec);
 684     AliTPCCalROC *oROC=fPulserOutlier->GetCalROC(isec);
 685     Float_t pt_mean=ptROC->GetMean(oROC);
 686     UInt_t nrows=mROC->GetNrows();
 687     for (UInt_t irow=0;irow<nrows;++irow){
 688       UInt_t npads=mROC->GetNPads(irow);
 689       for (UInt_t ipad=0;ipad<npads;++ipad){
 690         //don't use masked channels;
 691         if (mROC   ->GetValue(irow,ipad)) continue;
 692         if (mRefROC->GetValue(irow,ipad)) continue;
 693         //don't user edge pads
 694         if (ipad==0||ipad==npads-1) continue;
 695         //data
 696         Float_t pq=pqROC->GetValue(irow,ipad);
 697         Float_t pqRef=pqRefROC->GetValue(irow,ipad);
 698         Float_t pt=ptROC->GetValue(irow,ipad);
 699 //         Float_t ptRef=ptRefROC->GetValue(irow,ipad);
 700         //comparisons q
 701         Float_t deviation=TMath::Abs(pq/pqRef-1);
 702         for (Int_t i=0;i<npar;++i){
 703           if (deviation>vThres[i])
 704             ++pulserQdeviations.GetMatrixArray()[i];
 705         }
 706         if (pqRef>11&&pq<11) ++npadsOffAdd;
 707         varQMean+=pq-pqRef;
 708         //comparisons t
 709         if (TMath::Abs(pt-pt_mean)>1) ++npadsOutOneTB;
 710         ++nActive;
 711       }//end ipad
 712     }//ind irow
 713   }//end isec
 714   if (nActive>0){
 715     for (Int_t i=0;i<npar;++i){
 716       pulserQdeviations.GetMatrixArray()[i]/=nActive;
 717       varQMean/=nActive;
 718     }
 719   }
 720 }
 721 //_____________________________________________________________________________________
 722 void AliTPCcalibDButil::UpdatePulserOutlierMap()
 723 {
 724   //
 725   //
 726   //
 727   PulserOutlierMap(fPulserOutlier,fPulserTmean, fPulserQmean);
 728 }
 729 //_____________________________________________________________________________________
 730 void AliTPCcalibDButil::UpdateRefPulserOutlierMap()
 731 {
 732   //
 733   //
 734   //
 735   PulserOutlierMap(fRefPulserOutlier,fRefPulserTmean, fRefPulserQmean);
 736 }
 737 //_____________________________________________________________________________________
 738 void AliTPCcalibDButil::PulserOutlierMap(AliTPCCalPad *pulOut, const AliTPCCalPad *pulT, const AliTPCCalPad *pulQ)
 739 {
 740   //
 741   // Create a map that contains outliers from the Pulser calibration data.
 742   // The outliers include masked channels, edge pads and pads with
 743   //   too large timing and charge variations.
 744   // fNpulserOutliers is the number of outliers in the Pulser calibration data.
 745   //   those do not contain masked and edge pads
 746   //
 747   if (!pulT||!pulQ) {
 748     //reset map
 749     pulOut->Multiply(0.);
 750     fNpulserOutliers=-1;
 751     return;
 752   }
 753   AliTPCCalROC *rocMasked=0x0;
 754   fNpulserOutliers=0;
 755
 756   //Create Outlier Map
 757   for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 758     AliTPCCalROC *tmeanROC=pulT->GetCalROC(isec);
 759     AliTPCCalROC *qmeanROC=pulQ->GetCalROC(isec);
 760     AliTPCCalROC *outROC=pulOut->GetCalROC(isec);
 761     if (!tmeanROC||!qmeanROC) {
 762       //reset outliers in this ROC
 763       outROC->Multiply(0.);
 764       continue;
 765     }
 766     if (fALTROMasked) rocMasked=fALTROMasked->GetCalROC(isec);
 767 //     Double_t dummy=0;
 768 //     Float_t qmedian=qmeanROC->GetLTM(&dummy,.5);
 769 //     Float_t tmedian=tmeanROC->GetLTM(&dummy,.5);
 770     UInt_t nrows=tmeanROC->GetNrows();
 771     for (UInt_t irow=0;irow<nrows;++irow){
 772       UInt_t npads=tmeanROC->GetNPads(irow);
 773       for (UInt_t ipad=0;ipad<npads;++ipad){
 774         Int_t outlier=0,masked=0;
 775         Float_t q=qmeanROC->GetValue(irow,ipad);
 776         Float_t t=tmeanROC->GetValue(irow,ipad);
 777         //masked channels are outliers
 778         if (rocMasked && rocMasked->GetValue(irow,ipad)) masked=1;
 779         //edge pads are outliers
 780         if (ipad==0||ipad==npads-1) masked=1;
 781         //channels with too large charge or timing deviation from the meadian are outliers
 782 //         if (TMath::Abs(q-qmedian)>fPulQmaxLimitAbs || TMath::Abs(t-tmedian)>fPulTmaxLimitAbs) outlier=1;
 783         if (q<fPulQminLimit && !masked) outlier=1;
 784         //check for nan
 785         if ( !(q<10000000) || !(t<10000000)) outlier=1;
 786         outROC->SetValue(irow,ipad,outlier+masked);
 787         fNpulserOutliers+=outlier;
 788       }
 789     }
 790   }
 791 }
 792 //_____________________________________________________________________________________
 793 AliTPCCalPad* AliTPCcalibDButil::CreatePadTime0(Int_t model)
 794 {
 795   //
 796   // Create pad time0 object from pulser and/or CE data, depending on the selected model
 797   // Model 0: normalise each readout chamber to its mean, outlier cutted, only Pulser
 798   // Model 1: normalise IROCs/OROCs of each readout side to its mean, only Pulser
 799   //
 800   //
 801
 802   AliTPCCalPad *padTime0=new AliTPCCalPad("PadTime0",Form("PadTime0-Model_%d",model));
 803   // decide between different models
 804   if (model==0||model==1){
 805     TVectorD vMean;
 806     if (model==1) ProcessPulser(vMean);
 807     for (UInt_t isec=0;isec<AliTPCCalPad::kNsec;++isec){
 808       AliTPCCalROC *rocPulTmean=fPulserTmean->GetCalROC(isec);
 809       if (!rocPulTmean) continue;
 810       AliTPCCalROC *rocTime0=padTime0->GetCalROC(isec);
 811       AliTPCCalROC *rocOut=fPulserOutlier->GetCalROC(isec);
 812       Float_t mean=rocPulTmean->GetMean(rocOut);
 813       //treat case where a whole partition is masked
 814       if (mean==0) mean=rocPulTmean->GetMean();
 815       if (model==1) {
 816         Int_t type=isec/18;
 817         mean=vMean[type];
 818       }
 819       UInt_t nrows=rocTime0->GetNrows();
 820       for (UInt_t irow=0;irow<nrows;++irow){
 821         UInt_t npads=rocTime0->GetNPads(irow);
 822         for (UInt_t ipad=0;ipad<npads;++ipad){
 823           Float_t time=rocPulTmean->GetValue(irow,ipad);
 824           //in case of an outlier pad use the mean of the altro values.
 825           //This should be the most precise guess in that case.
 826           if (rocOut->GetValue(irow,ipad)) {
 827             time=GetMeanAltro(rocPulTmean,irow,ipad,rocOut);
 828             if (time==0) time=mean;
 829           }
 830           Float_t val=time-mean;
 831           rocTime0->SetValue(irow,ipad,val);
 832         }
 833       }
 834     }
 835   }
 836
 837
 838
 839   return padTime0;
 840 }
 841 //_____________________________________________________________________________________
 842 Float_t AliTPCcalibDButil::GetMeanAltro(const AliTPCCalROC *roc, const Int_t row, const Int_t pad, AliTPCCalROC *rocOut)
 843 {
 844   if (roc==0) return 0.;
 845   const Int_t sector=roc->GetSector();
 846   AliTPCROC *tpcRoc=AliTPCROC::Instance();
 847   const UInt_t altroRoc=fMapper->GetFEC(sector,row,pad)*8+fMapper->GetChip(sector,row,pad);
 848   Float_t mean=0;
 849   Int_t   n=0;
 850
 851   //loop over a small range around the requested pad (+-10 rows/pads)
 852   for (Int_t irow=row-10;irow<row+10;++irow){
 853     if (irow<0||irow>(Int_t)tpcRoc->GetNRows(sector)-1) continue;
 854     for (Int_t ipad=pad-10; ipad<pad+10;++ipad){
 855       if (ipad<0||ipad>(Int_t)tpcRoc->GetNPads(sector,irow)-1) continue;
 856       const UInt_t altroCurr=fMapper->GetFEC(sector,irow,ipad)*8+fMapper->GetChip(sector,irow,ipad);
 857       if (altroRoc!=altroCurr) continue;
 858       if ( rocOut && rocOut->GetValue(irow,ipad) ) continue;
 859       Float_t val=roc->GetValue(irow,ipad);
 860       mean+=val;
 861       ++n;
 862     }
 863   }
 864   if (n>0) mean/=n;
 865   return mean;
 866 }
 867 //_____________________________________________________________________________________
 868 void AliTPCcalibDButil::SetRefFile(const char* filename)
 869 {
 870   //
 871   // load cal pad objects form the reference file
 872   //
 873   TDirectory *currDir=gDirectory;
 874   TFile f(filename);
 875   fRefPedestals=(AliTPCCalPad*)f.Get("Pedestals");
 876   fRefPadNoise=(AliTPCCalPad*)f.Get("PadNoise");
 877   //pulser data
 878   fRefPulserTmean=(AliTPCCalPad*)f.Get("PulserTmean");
 879   fRefPulserTrms=(AliTPCCalPad*)f.Get("PulserTrms");
 880   fRefPulserQmean=(AliTPCCalPad*)f.Get("PulserQmean");
 881   //CE data
 882   fRefCETmean=(AliTPCCalPad*)f.Get("CETmean");
 883   fRefCETrms=(AliTPCCalPad*)f.Get("CETrms");
 884   fRefCEQmean=(AliTPCCalPad*)f.Get("CEQmean");
 885   //Altro data
 886 //   fRefALTROAcqStart=(AliTPCCalPad*)f.Get("ALTROAcqStart");
 887 //   fRefALTROZsThr=(AliTPCCalPad*)f.Get("ALTROZsThr");
 888 //   fRefALTROFPED=(AliTPCCalPad*)f.Get("ALTROFPED");
 889 //   fRefALTROAcqStop=(AliTPCCalPad*)f.Get("ALTROAcqStop");
 890   fRefALTROMasked=(AliTPCCalPad*)f.Get("ALTROMasked");
 891   f.Close();
 892   currDir->cd();
 893 }
 894