TPC/AliTPCCalibPulser.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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////////////////////////////////
  19 //                                                                                     //
  20 //                  Implementation of the TPC pulser calibration                       //
  21 //                                                                                     //
  22 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch   //
  23 //                                                                                     //
  24 /////////////////////////////////////////////////////////////////////////////////////////
  25 /***************************************************************************
  26  *                      Class Description                                  *
  27  ***************************************************************************
  28
  29  The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
  30  runs performed with the calibration pulser.
  31
  32  The information retrieved is
  33  - Time0 differences
  34  - Signal width differences
  35  - Amplification variations
  36
  37  the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
  38  one chip and somewhat large between different chips.
  39
  40  Histograms:
  41    For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum]) is created when
  42    it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
  43    TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
  44
  45
  46  Working principle:
  47  ------------------
  48  Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
  49  (see below). These in the end call the Update(...) function.
  50
  51  - the Update(...) function:
  52    In this function the array fPadSignal is filled with the adc signals between the specified range
  53    fFirstTimeBin and fLastTimeBin for the current pad.
  54    before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
  55    stored in fPadSignal.
  56
  57    - the ProcessPad() function:
  58      Find Pedestal and Noise information
  59      - use database information which has to be set by calling
  60        SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
  61      - if no information from the pedestal data base
  62        is available the informaion is calculated on the fly ( see FindPedestal() function )
  63
  64      Find the Pulser signal information
  65      - calculate  mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
  66        the Q sum is scaled by pad area
  67        (see FindPulserSignal(...) function)
  68
  69      Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
  70      Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
  71
  72  At the end of each event the EndEvent() function is called
  73
  74  - the EndEvent() function:
  75    calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
  76    This is done to overcome syncronisation problems between the trigger and the fec clock.
  77
  78  After accumulating the desired statistics the Analyse() function has to be called.
  79  - the Analyse() function
  80    Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
  81    the AliMathBase::GetCOG(...) function, and the calibration
  82    storage classes (AliTPCCalROC) are filled for each ROC.
  83    The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
  84    fCalRocArrayQ;
  85
  86
  87
  88  User interface for filling data:
  89  --------------------------------
  90
  91  To Fill information one of the following functions can be used:
  92
  93  Bool_t ProcessEvent(eventHeaderStruct *event);
  94    - process Date event
  95    - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
  96
  97  Bool_t ProcessEvent(AliRawReader *rawReader);
  98    - process AliRawReader event
  99    - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
 100
 101  Bool_t ProcessEvent(AliTPCRawStream *rawStream);
 102    - process event from AliTPCRawStream
 103    - call Update function for signal filling
 104
 105  Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 106               iPad,  const Int_t iTimeBin, const Float_t signal);
 107    - directly  fill signal information (sector, row, pad, time bin, pad)
 108      to the reference histograms
 109
 110  It is also possible to merge two independently taken calibrations using the function
 111
 112  void Merge(AliTPCCalibPulser *sig)
 113    - copy histograms in 'sig' if the do not exist in this instance
 114    - Add histograms in 'sig' to the histograms in this instance if the allready exist
 115    - After merging call Analyse again!
 116
 117
 118
 119  -- example: filling data using root raw data:
 120  void fillSignal(Char_t *filename)
 121  {
 122     rawReader = new AliRawReaderRoot(fileName);
 123     if ( !rawReader ) return;
 124     AliTPCCalibPulser *calib = new AliTPCCalibPulser;
 125     while (rawReader->NextEvent()){
 126       calib->ProcessEvent(rawReader);
 127     }
 128     calib->Analyse();
 129     calib->DumpToFile("SignalData.root");
 130     delete rawReader;
 131     delete calib;
 132  }
 133
 134
 135  What kind of information is stored and how to retrieve them:
 136  ------------------------------------------------------------
 137
 138  - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
 139
 140    TH2F *GetHistoT0(Int_t sector);
 141    TH2F *GetHistoRMS(Int_t sector);
 142    TH2F *GetHistoQ(Int_t sector);
 143
 144  - Accessing the calibration storage objects:
 145
 146    AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values
 147    AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values
 148    AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values
 149
 150    example for visualisation:
 151    if the file "SignalData.root" was created using the above example one could do the following:
 152
 153    TFile fileSignal("SignalData.root")
 154    AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
 155    sig->GetCalRocT0(0)->Draw("colz");
 156    sig->GetCalRocRMS(0)->Draw("colz");
 157
 158    or use the AliTPCCalPad functionality:
 159    AliTPCCalPad padT0(ped->GetCalPadT0());
 160    AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
 161    padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
 162    padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
 163 */
 164
 165
 166 //Root includes
 167 #include <TObjArray.h>
 168 #include <TH1F.h>
 169 #include <TH2F.h>
 170 #include <TH2S.h>
 171 #include <TString.h>
 172 #include <TVectorF.h>
 173 #include <TMath.h>
 174 #include <TMap.h>
 175 #include <TDirectory.h>
 176 #include <TSystem.h>
 177 #include <TROOT.h>
 178 #include <TFile.h>
 179
 180 //AliRoot includes
 181 #include "AliRawReader.h"
 182 #include "AliRawReaderRoot.h"
 183 #include "AliRawReaderDate.h"
 184 #include "AliTPCRawStream.h"
 185 #include "AliTPCRawStreamFast.h"
 186 #include "AliTPCCalROC.h"
 187 #include "AliTPCCalPad.h"
 188 #include "AliTPCROC.h"
 189 #include "AliTPCParam.h"
 190 #include "AliTPCCalibPulser.h"
 191 #include "AliTPCcalibDB.h"
 192 #include "AliMathBase.h"
 193 #include "AliLog.h"
 194 #include "TTreeStream.h"
 195
 196 //date
 197 #include "event.h"
 198
 199
 200
 201
 202 ClassImp(AliTPCCalibPulser)
 203
 204 AliTPCCalibPulser::AliTPCCalibPulser() :
 205 AliTPCCalibRawBase(),
 206 fNbinsT0(200),
 207 fXminT0(-2),
 208 fXmaxT0(2),
 209 fNbinsQ(200),
 210 fXminQ(10),
 211 fXmaxQ(40),
 212 fNbinsRMS(100),
 213 fXminRMS(0.1),
 214 fXmaxRMS(5.1),
 215 fPeakIntMinus(2),
 216 fPeakIntPlus(2),
 217 fIsZeroSuppressed(kFALSE),
 218 fLastSector(-1),
 219 fParam(new AliTPCParam),
 220 fPedestalTPC(0x0),
 221 fPadNoiseTPC(0x0),
 222 fOutliers(0x0),
 223 fPedestalROC(0x0),
 224 fPadNoiseROC(0x0),
 225 fCalRocArrayT0(72),
 226 fCalRocArrayQ(72),
 227 fCalRocArrayRMS(72),
 228 fCalRocArrayOutliers(72),
 229 fHistoQArray(72),
 230 fHistoT0Array(72),
 231 fHistoRMSArray(72),
 232 fHMeanTimeSector(0x0),
 233 fVMeanTimeSector(72),
 234 fPadTimesArrayEvent(72),
 235 fPadQArrayEvent(72),
 236 fPadRMSArrayEvent(72),
 237 fPadPedestalArrayEvent(72),
 238 fCurrentChannel(-1),
 239 fCurrentSector(-1),
 240 fCurrentRow(-1),
 241 fCurrentPad(-1),
 242 fMaxPadSignal(-1),
 243 fMaxTimeBin(-1),
 244 fPadSignal(1024),
 245 fPadPedestal(0),
 246 fPadNoise(0),
 247 fVTime0Offset(72),
 248 fVTime0OffsetCounter(72)
 249 {
 250   //
 251   // AliTPCSignal default constructor
 252   //
 253   SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
 254   fFirstTimeBin=60;
 255   fLastTimeBin=900;
 256   fParam->Update();
 257 }
 258 //_____________________________________________________________________
 259 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
 260 AliTPCCalibRawBase(sig),
 261 fNbinsT0(sig.fNbinsT0),
 262 fXminT0(sig.fXminT0),
 263 fXmaxT0(sig.fXmaxT0),
 264 fNbinsQ(sig.fNbinsQ),
 265 fXminQ(sig.fXminQ),
 266 fXmaxQ(sig.fXmaxQ),
 267 fNbinsRMS(sig.fNbinsRMS),
 268 fXminRMS(sig.fXminRMS),
 269 fXmaxRMS(sig.fXmaxRMS),
 270 fPeakIntMinus(sig.fPeakIntMinus),
 271 fPeakIntPlus(sig.fPeakIntPlus),
 272 fIsZeroSuppressed(sig.fIsZeroSuppressed),
 273 fLastSector(-1),
 274 fParam(new AliTPCParam),
 275 fPedestalTPC(0x0),
 276 fPadNoiseTPC(0x0),
 277 fOutliers(0x0),
 278 fPedestalROC(0x0),
 279 fPadNoiseROC(0x0),
 280 fCalRocArrayT0(72),
 281 fCalRocArrayQ(72),
 282 fCalRocArrayRMS(72),
 283 fCalRocArrayOutliers(72),
 284 fHistoQArray(72),
 285 fHistoT0Array(72),
 286 fHistoRMSArray(72),
 287 fHMeanTimeSector(0x0),
 288 fVMeanTimeSector(72),
 289 fPadTimesArrayEvent(72),
 290 fPadQArrayEvent(72),
 291 fPadRMSArrayEvent(72),
 292 fPadPedestalArrayEvent(72),
 293 fCurrentChannel(-1),
 294 fCurrentSector(-1),
 295 fCurrentRow(-1),
 296 fCurrentPad(-1),
 297 fMaxPadSignal(-1),
 298 fMaxTimeBin(-1),
 299 fPadSignal(1024),
 300 fPadPedestal(0),
 301 fPadNoise(0),
 302 fVTime0Offset(72),
 303 fVTime0OffsetCounter(72)
 304 {
 305   //
 306   // AliTPCSignal default constructor
 307   //
 308
 309   for (Int_t iSec = 0; iSec < 72; ++iSec){
 310     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 311     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 312     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 313     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 314
 315     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 316     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 317     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 318
 319     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 320     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 321     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 322     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 323
 324     if ( hQ != 0x0 ){
 325       TH2S *hNew = new TH2S(*hQ);
 326       hNew->SetDirectory(0);
 327       fHistoQArray.AddAt(hNew,iSec);
 328     }
 329     if ( hT0 != 0x0 ){
 330       TH2S *hNew = new TH2S(*hT0);
 331       hNew->SetDirectory(0);
 332       fHistoQArray.AddAt(hNew,iSec);
 333     }
 334     if ( hRMS != 0x0 ){
 335       TH2S *hNew = new TH2S(*hRMS);
 336       hNew->SetDirectory(0);
 337       fHistoQArray.AddAt(hNew,iSec);
 338     }
 339     fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
 340   }
 341
 342   if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
 343   fParam->Update();
 344 }
 345 AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
 346 AliTPCCalibRawBase(),
 347 fNbinsT0(200),
 348 fXminT0(-2),
 349 fXmaxT0(2),
 350 fNbinsQ(200),
 351 fXminQ(10),
 352 fXmaxQ(40),
 353 fNbinsRMS(100),
 354 fXminRMS(0.1),
 355 fXmaxRMS(5.1),
 356 fPeakIntMinus(2),
 357 fPeakIntPlus(2),
 358 fIsZeroSuppressed(kFALSE),
 359 fLastSector(-1),
 360 fParam(new  AliTPCParam),
 361 fPedestalTPC(0x0),
 362 fPadNoiseTPC(0x0),
 363 fOutliers(0x0),
 364 fPedestalROC(0x0),
 365 fPadNoiseROC(0x0),
 366 fCalRocArrayT0(72),
 367 fCalRocArrayQ(72),
 368 fCalRocArrayRMS(72),
 369 fCalRocArrayOutliers(72),
 370 fHistoQArray(72),
 371 fHistoT0Array(72),
 372 fHistoRMSArray(72),
 373 fHMeanTimeSector(0x0),
 374 fVMeanTimeSector(72),
 375 fPadTimesArrayEvent(72),
 376 fPadQArrayEvent(72),
 377 fPadRMSArrayEvent(72),
 378 fPadPedestalArrayEvent(72),
 379 fCurrentChannel(-1),
 380 fCurrentSector(-1),
 381 fCurrentRow(-1),
 382 fCurrentPad(-1),
 383 fMaxPadSignal(-1),
 384 fMaxTimeBin(-1),
 385 fPadSignal(1024),
 386 fPadPedestal(0),
 387 fPadNoise(0),
 388 fVTime0Offset(72),
 389 fVTime0OffsetCounter(72)
 390 {
 391   //
 392   // This constructor uses a TMap for setting some parametes
 393   //
 394   SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
 395   fFirstTimeBin=60;
 396   fLastTimeBin=900;
 397   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 398   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 399   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 400   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 401   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 402   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 403   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 404   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 405   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 406   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 407   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 408   if (config->GetValue("PeakIntMinus")) fPeakIntMinus = (Int_t)((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atof();
 409   if (config->GetValue("PeakIntPlus")) fPeakIntPlus = (Int_t)((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atof();
 410   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 411
 412   fParam->Update();
 413 }
 414 //_____________________________________________________________________
 415 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
 416 {
 417   //
 418   // assignment operator
 419   //
 420   if (&source == this) return *this;
 421   new (this) AliTPCCalibPulser(source);
 422
 423   return *this;
 424 }
 425 //_____________________________________________________________________
 426 AliTPCCalibPulser::~AliTPCCalibPulser()
 427 {
 428   //
 429   // destructor
 430   //
 431
 432   Reset();
 433   delete fParam;
 434 }
 435 void AliTPCCalibPulser::Reset()
 436 {
 437   //
 438   // Delete all information: Arrays, Histograms, CalRoc objects
 439   //
 440   fCalRocArrayT0.Delete();
 441   fCalRocArrayQ.Delete();
 442   fCalRocArrayRMS.Delete();
 443   fCalRocArrayOutliers.Delete();
 444
 445   fHistoQArray.Delete();
 446   fHistoT0Array.Delete();
 447   fHistoRMSArray.Delete();
 448
 449   fPadTimesArrayEvent.Delete();
 450   fPadQArrayEvent.Delete();
 451   fPadRMSArrayEvent.Delete();
 452   fPadPedestalArrayEvent.Delete();
 453
 454   if (fHMeanTimeSector) delete fHMeanTimeSector;
 455 }
 456 //_____________________________________________________________________
 457 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
 458                                 const Int_t icRow,
 459                                 const Int_t icPad,
 460                                 const Int_t icTimeBin,
 461                                 const Float_t csignal)
 462 {
 463     //
 464     // Signal filling methode on the fly pedestal and time offset correction if necessary.
 465     // no extra analysis necessary. Assumes knowledge of the signal shape!
 466     // assumes that it is looped over consecutive time bins of one pad
 467     //
 468
 469   if (icRow<0) return 0;
 470   if (icPad<0) return 0;
 471   if (icTimeBin<0) return 0;
 472   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 473
 474   if ( icRow<0 || icPad<0 ){
 475     AliWarning("Wrong Pad or Row number, skipping!");
 476     return 0;
 477   }
 478
 479   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 480
 481     //init first pad and sector in this event
 482   if ( fCurrentChannel == -1 ) {
 483     fCurrentChannel = iChannel;
 484     fCurrentSector  = icsector;
 485     fCurrentRow     = icRow;
 486     fCurrentPad     = icPad;
 487   }
 488
 489     //process last pad if we change to a new one
 490   if ( iChannel != fCurrentChannel ){
 491     ProcessPad();
 492     fLastSector=fCurrentSector;
 493     fCurrentChannel = iChannel;
 494     fCurrentSector  = icsector;
 495     fCurrentRow     = icRow;
 496     fCurrentPad     = icPad;
 497   }
 498
 499     //fill signals for current pad
 500   fPadSignal[icTimeBin]=csignal;
 501   if ( csignal > fMaxPadSignal ){
 502     fMaxPadSignal = csignal;
 503     fMaxTimeBin   = icTimeBin;
 504   }
 505   return 0;
 506 }
 507 //_____________________________________________________________________
 508 void AliTPCCalibPulser::FindPedestal(Float_t part)
 509 {
 510     //
 511     // find pedestal and noise for the current pad. Use either database or
 512     // truncated mean with part*100%
 513     //
 514   Bool_t noPedestal = kTRUE;;
 515   if (fPedestalTPC&&fPadNoiseTPC){
 516         //use pedestal database
 517         //only load new pedestals if the sector has changed
 518     if ( fCurrentSector!=fLastSector ){
 519       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 520       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 521     }
 522
 523     if ( fPedestalROC&&fPadNoiseROC ){
 524       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 525       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 526       noPedestal   = kFALSE;
 527     }
 528
 529   }
 530
 531     //if we are not running with pedestal database, or for the current sector there is no information
 532     //available, calculate the pedestal and noise on the fly
 533   if ( noPedestal ) {
 534     const Int_t kPedMax = 100;  //maximum pedestal value
 535     Float_t  max    =  0;
 536     Float_t  maxPos =  0;
 537     Int_t    median =  -1;
 538     Int_t    count0 =  0;
 539     Int_t    count1 =  0;
 540   //
 541     Float_t padSignal=0;
 542         //
 543     UShort_t histo[kPedMax];
 544     memset(histo,0,kPedMax*sizeof(UShort_t));
 545
 546     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 547       padSignal = fPadSignal.GetMatrixArray()[i];
 548       if (padSignal<=0) continue;
 549       if (padSignal>max && i>10) {
 550         max = padSignal;
 551         maxPos = i;
 552       }
 553       if (padSignal>kPedMax-1) continue;
 554       histo[Int_t(padSignal+0.5)]++;
 555       count0++;
 556     }
 557       //
 558     for (Int_t i=1; i<kPedMax; ++i){
 559       if (count1<count0*0.5) median=i;
 560       count1+=histo[i];
 561     }
 562   // truncated mean
 563   //
 564         // what if by chance histo[median] == 0 ?!?
 565     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 566   //
 567     for (Int_t idelta=1; idelta<10; ++idelta){
 568       if (median-idelta<=0) continue;
 569       if (median+idelta>kPedMax) continue;
 570       if (count<part*count1){
 571         count+=histo[median-idelta];
 572         mean +=histo[median-idelta]*(median-idelta);
 573         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 574         count+=histo[median+idelta];
 575         mean +=histo[median+idelta]*(median+idelta);
 576         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 577       }
 578     }
 579     fPadPedestal = 0;
 580     fPadNoise    = 0;
 581     if ( count > 0 ) {
 582       mean/=count;
 583       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 584       fPadPedestal = mean;
 585       fPadNoise    = rms;
 586     }
 587   }
 588   fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
 589 }
 590 //_____________________________________________________________________
 591 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
 592 {
 593 //
 594     //  Find position, signal width and height of the CE signal (last signal)
 595     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 596     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 597     //
 598
 599   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 600   Int_t   cemaxpos       = fMaxTimeBin;
 601   Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal sum
 602   Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal max
 603   const Int_t    kCemin  = fPeakIntMinus;             // range for the analysis of the ce signal +- channels from the peak
 604   const Int_t    kCemax  = fPeakIntPlus;
 605   param[0] = ceQmax;
 606   param[1] = ceTime;
 607   param[2] = ceRMS;
 608   qSum     = ceQsum;
 609
 610   if (cemaxpos>0){
 611     ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 612     if ( ceQmax<ceMaxThreshold ) return;
 613     for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
 614       Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 615       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 616         ceTime+=signal*(i+0.5);
 617         ceRMS +=signal*(i+0.5)*(i+0.5);
 618         ceQsum+=signal;
 619       }
 620     }
 621   }
 622   if (ceQsum>ceSumThreshold) {
 623     ceTime/=ceQsum;
 624     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 625     ceTime-=GetL1PhaseTB();
 626         //only fill the Time0Offset if pad was not marked as an outlier!
 627     if ( !fOutliers ){
 628       //skip edge pads for calculating the mean time
 629       if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
 630         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 631         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 632         GetHistoTSec()->Fill(ceTime,fCurrentSector);
 633       }
 634     } else {
 635       if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
 636         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 637         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 638       }
 639     }
 640
 641   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
 642   //                                the pick-up signal should scale with the pad area. In addition
 643   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
 644   //                                ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
 645   //                                to the OROC signal (factor 2/3 for the IROCs).
 646     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
 647     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
 648
 649     ceQsum/=norm;
 650   } else {
 651     ceQmax=0;
 652     ceTime=0;
 653     ceRMS =0;
 654     ceQsum=0;
 655   }
 656   param[0] = ceQmax;
 657   param[1] = ceTime;
 658   param[2] = ceRMS;
 659   qSum     = ceQsum;
 660 }
 661 //_____________________________________________________________________
 662 void AliTPCCalibPulser::ProcessPad()
 663 {
 664   //
 665   //  Process data of current pad
 666   //
 667
 668   FindPedestal();
 669   TVectorD param(3);
 670   Float_t  qSum;
 671   FindPulserSignal(param, qSum);
 672
 673   Double_t meanT  = param[1];
 674   Double_t sigmaT = param[2];
 675
 676
 677   //Fill Event T0 counter
 678   (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 679
 680   //Fill Q histogram
 681 //   GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
 682   GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
 683
 684   //Fill RMS histogram
 685   GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 686
 687
 688     //Fill debugging info
 689   if ( GetStreamLevel()>0 ){
 690     TTreeSRedirector *streamer=GetDebugStreamer();
 691     if ( GetStreamLevel() == 1 ){
 692       if ( streamer ) {
 693         Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
 694         (*streamer) << "PadSignals" <<
 695           "Event="  <<fNevents <<
 696           "Sector=" <<fCurrentSector<<
 697           "Row="    <<fCurrentRow<<
 698           "Pad="    <<fCurrentPad<<
 699           "PadC="   <<padc<<
 700           "Channel="<<fCurrentChannel<<
 701           "Sum="    <<qSum<<
 702           "params.="<<&param<<
 703           "signal.=" <<&fPadSignal<<
 704           "\n";
 705       }
 706     } else { //debug > 1
 707       (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
 708       (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 709       (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
 710     }
 711   }
 712   ResetPad();
 713 }
 714 //_____________________________________________________________________
 715 void AliTPCCalibPulser::EndEvent()
 716 {
 717     //
 718     //  Process data of current event
 719     //
 720
 721     //check if last pad has allready been processed, if not do so
 722   if ( fMaxTimeBin>-1 ) ProcessPad();
 723
 724     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 725   for ( Int_t iSec = 0; iSec<72; ++iSec ){
 726     TVectorF *vTimes = GetPadTimesEvent(iSec);
 727     if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
 728     Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 729     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 730       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
 731
 732       GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
 733             //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
 734
 735
 736       //Debug start
 737       if ( GetStreamLevel()>1 ){
 738         TTreeSRedirector *streamer=GetDebugStreamer();
 739         if ( streamer ) {
 740           Int_t row=0;
 741           Int_t pad=0;
 742           Int_t padc=0;
 743
 744           Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 745           Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 746
 747           UInt_t channel=iChannel;
 748           Int_t sector=iSec;
 749
 750           while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 751           pad = channel-fROC->GetRowIndexes(sector)[row];
 752           padc = pad-(fROC->GetNPads(sector,row)/2);
 753
 754           (*streamer) << "DataPad" <<
 755                     "Event=" << fNevents <<
 756             "Sector="<< sector <<
 757             "Row="   << row<<
 758             "Pad="   << pad <<
 759             "PadC="  << padc <<
 760             "PadSec="<< channel <<
 761             "Time0="  << time0 <<
 762             "Time="  << time <<
 763             "RMS="   << rms <<
 764             "Sum="   << q <<
 765             "\n";
 766         }
 767       }
 768       //Debug end
 769     }
 770   }
 771   ++fNevents;
 772 }
 773 //_____________________________________________________________________
 774 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
 775                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 776                                   const Char_t *type, Bool_t force)
 777 {
 778     //
 779     // return pointer to Q histogram
 780     // if force is true create a new histogram if it doesn't exist allready
 781     //
 782   if ( !force || arr->UncheckedAt(sector) )
 783     return (TH2S*)arr->UncheckedAt(sector);
 784
 785   // if we are forced and histogram doesn't yes exist create it
 786   // new histogram with Q calib information. One value for each pad!
 787   TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
 788                         nbinsY, ymin, ymax,
 789                         fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 790   hist->SetDirectory(0);
 791   arr->AddAt(hist,sector);
 792   return hist;
 793 }
 794 //_____________________________________________________________________
 795 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
 796 {
 797     //
 798     // return pointer to T0 histogram
 799     // if force is true create a new histogram if it doesn't exist allready
 800     //
 801   TObjArray *arr = &fHistoT0Array;
 802   return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 803 }
 804 //_____________________________________________________________________
 805 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
 806 {
 807     //
 808     // return pointer to Q histogram
 809     // if force is true create a new histogram if it doesn't exist allready
 810     //
 811   TObjArray *arr = &fHistoQArray;
 812   return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 813 }
 814 //_____________________________________________________________________
 815 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
 816 {
 817     //
 818     // return pointer to Q histogram
 819     // if force is true create a new histogram if it doesn't exist allready
 820     //
 821   TObjArray *arr = &fHistoRMSArray;
 822   return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 823 }
 824 //_____________________________________________________________________
 825 TH2F* AliTPCCalibPulser::GetHistoTSec()
 826 {
 827     //
 828     // return the pointer to the abs time distribution per sector
 829     // create it if it does not exist
 830     //
 831   if ( !fHMeanTimeSector )   //!!!if you change the binning here, you should also change it in the Analyse Function!!
 832     fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
 833                                 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
 834                                 72,0,72);
 835   return fHMeanTimeSector;
 836 }
 837 //_____________________________________________________________________
 838 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
 839 {
 840     //
 841     // return pointer to Pad Info from 'arr' for the current event and sector
 842     // if force is true create it if it doesn't exist allready
 843     //
 844   if ( !force || arr->UncheckedAt(sector) )
 845     return (TVectorF*)arr->UncheckedAt(sector);
 846
 847   TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 848   arr->AddAt(vect,sector);
 849   return vect;
 850 }
 851 //_____________________________________________________________________
 852 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
 853 {
 854     //
 855     // return pointer to Pad Times Array for the current event and sector
 856     // if force is true create it if it doesn't exist allready
 857     //
 858   TObjArray *arr = &fPadTimesArrayEvent;
 859   return GetPadInfoEvent(sector,arr,force);
 860 }
 861 //_____________________________________________________________________
 862 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
 863 {
 864     //
 865     // return pointer to Pad Q Array for the current event and sector
 866     // if force is true create it if it doesn't exist allready
 867     // for debugging purposes only
 868     //
 869
 870   TObjArray *arr = &fPadQArrayEvent;
 871   return GetPadInfoEvent(sector,arr,force);
 872 }
 873 //_____________________________________________________________________
 874 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
 875 {
 876     //
 877     // return pointer to Pad RMS Array for the current event and sector
 878     // if force is true create it if it doesn't exist allready
 879     // for debugging purposes only
 880     //
 881   TObjArray *arr = &fPadRMSArrayEvent;
 882   return GetPadInfoEvent(sector,arr,force);
 883 }
 884 //_____________________________________________________________________
 885 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
 886 {
 887     //
 888     // return pointer to Pad RMS Array for the current event and sector
 889     // if force is true create it if it doesn't exist allready
 890     // for debugging purposes only
 891     //
 892   TObjArray *arr = &fPadPedestalArrayEvent;
 893   return GetPadInfoEvent(sector,arr,force);
 894 }
 895 //_____________________________________________________________________
 896 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
 897 {
 898     //
 899     // return pointer to ROC Calibration
 900     // if force is true create a new histogram if it doesn't exist allready
 901     //
 902   if ( !force || arr->UncheckedAt(sector) )
 903     return (AliTPCCalROC*)arr->UncheckedAt(sector);
 904
 905     // if we are forced and histogram doesn't yes exist create it
 906
 907     // new AliTPCCalROC for T0 information. One value for each pad!
 908   AliTPCCalROC *croc = new AliTPCCalROC(sector);
 909   arr->AddAt(croc,sector);
 910   return croc;
 911 }
 912 //_____________________________________________________________________
 913 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
 914 {
 915     //
 916     // return pointer to Carge ROC Calibration
 917     // if force is true create a new histogram if it doesn't exist allready
 918     //
 919   TObjArray *arr = &fCalRocArrayT0;
 920   return GetCalRoc(sector, arr, force);
 921 }
 922 //_____________________________________________________________________
 923 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
 924 {
 925     //
 926     // return pointer to T0 ROC Calibration
 927     // if force is true create a new histogram if it doesn't exist allready
 928     //
 929   TObjArray *arr = &fCalRocArrayQ;
 930   return GetCalRoc(sector, arr, force);
 931 }
 932 //_____________________________________________________________________
 933 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
 934 {
 935     //
 936     // return pointer to signal width ROC Calibration
 937     // if force is true create a new histogram if it doesn't exist allready
 938     //
 939   TObjArray *arr = &fCalRocArrayRMS;
 940   return GetCalRoc(sector, arr, force);
 941 }
 942 //_____________________________________________________________________
 943 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 944 {
 945     //
 946     // return pointer to Outliers
 947     // if force is true create a new histogram if it doesn't exist allready
 948     //
 949   TObjArray *arr = &fCalRocArrayOutliers;
 950   return GetCalRoc(sector, arr, force);
 951 }
 952 //_____________________________________________________________________
 953 void AliTPCCalibPulser::ResetEvent()
 954 {
 955     //
 956     //  Reset global counters  -- Should be called before each event is processed
 957     //
 958   fLastSector=-1;
 959   fCurrentSector=-1;
 960   fCurrentRow=-1;
 961   fCurrentPad=-1;
 962   fCurrentChannel=-1;
 963
 964   ResetPad();
 965
 966   fPadTimesArrayEvent.Delete();
 967
 968   fPadQArrayEvent.Delete();
 969   fPadRMSArrayEvent.Delete();
 970   fPadPedestalArrayEvent.Delete();
 971
 972   for ( Int_t i=0; i<72; ++i ){
 973     fVTime0Offset[i]=0;
 974     fVTime0OffsetCounter[i]=0;
 975   }
 976 }
 977 //_____________________________________________________________________
 978 void AliTPCCalibPulser::ResetPad()
 979 {
 980     //
 981     //  Reset pad infos -- Should be called after a pad has been processed
 982     //
 983   for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 984     fPadSignal[i] = 0;
 985   fMaxTimeBin = -1;
 986   fMaxPadSignal = -1;
 987   fPadPedestal  = -1;
 988   fPadNoise     = -1;
 989 }
 990 //_____________________________________________________________________
 991 Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
 992 {
 993     //
 994     // return true if pad is on the edge of a row
 995     //
 996   Int_t edge1   = 0;
 997   Int_t edge2   = fROC->GetNPads(sector,row)-1;
 998   if ( pad == edge1 || pad == edge2 ) return kTRUE;
 999
1000   return kFALSE;
1001 }
1002 //_____________________________________________________________________
1003 void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
1004 {
1005   //
1006   //  Merge reference histograms of sig to the current AliTPCCalibPulser
1007   //
1008
1009   MergeBase(sig);
1010   //merge histograms
1011   for (Int_t iSec=0; iSec<72; ++iSec){
1012     TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
1013     TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
1014     TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1015
1016
1017     if ( hRefQmerge ){
1018       TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1019       TH2S *hRefQ   = GetHistoQ(iSec);
1020       if ( hRefQ ) hRefQ->Add(hRefQmerge);
1021       else {
1022         TH2S *hist = new TH2S(*hRefQmerge);
1023         hist->SetDirectory(0);
1024         fHistoQArray.AddAt(hist, iSec);
1025       }
1026       hRefQmerge->SetDirectory(dir);
1027     }
1028     if ( hRefT0merge ){
1029       TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1030       TH2S *hRefT0  = GetHistoT0(iSec);
1031       if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1032       else {
1033         TH2S *hist = new TH2S(*hRefT0merge);
1034         hist->SetDirectory(0);
1035         fHistoT0Array.AddAt(hist, iSec);
1036       }
1037       hRefT0merge->SetDirectory(dir);
1038     }
1039     if ( hRefRMSmerge ){
1040       TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1041       TH2S *hRefRMS = GetHistoRMS(iSec);
1042       if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1043       else {
1044         TH2S *hist = new TH2S(*hRefRMSmerge);
1045         hist->SetDirectory(0);
1046         fHistoRMSArray.AddAt(hist, iSec);
1047       }
1048       hRefRMSmerge->SetDirectory(dir);
1049     }
1050
1051   }
1052   if ( sig->fHMeanTimeSector ){
1053     TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1054     if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1055     else {
1056       fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1057       fHMeanTimeSector->SetDirectory(0);
1058     }
1059     sig->fHMeanTimeSector->SetDirectory(dir);
1060   }
1061 }
1062
1063
1064 //_____________________________________________________________________
1065 Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1066 {
1067   //
1068   // Merge all objects of this type in list
1069   //
1070
1071   Long64_t nmerged=1;
1072
1073   TIter next(list);
1074   AliTPCCalibPulser *ce=0;
1075   TObject *o=0;
1076
1077   while ( (o=next()) ){
1078     ce=dynamic_cast<AliTPCCalibPulser*>(o);
1079     if (ce){
1080       Merge(ce);
1081       ++nmerged;
1082     }
1083   }
1084
1085   return nmerged;
1086 }
1087
1088 //_____________________________________________________________________
1089 void AliTPCCalibPulser::Analyse()
1090 {
1091   //
1092   //  Calculate calibration constants
1093   //
1094
1095   TVectorD paramQ(3);
1096   TVectorD paramT0(3);
1097   TVectorD paramRMS(3);
1098   TMatrixD dummy(3,3);
1099   //calculate mean time for each sector and mean time for each side
1100   TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1101   fVMeanTimeSector.Zero();
1102
1103   for (Int_t iSec=0; iSec<72; ++iSec){
1104     TH2S *hT0 = GetHistoT0(iSec);
1105     if (!hT0 ) continue;
1106     //calculate sector mean T
1107     if ( fHMeanTimeSector ){
1108       Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1109       Int_t offset = (nbinsT+2)*(iSec+1);
1110       Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1111       Int_t entries=0;
1112       for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1113       hMeanTsec.Set(nbinsT+2,arrP);
1114       hMeanTsec.SetEntries(entries);
1115       paramT0.Zero();
1116       // truncated mean: remove lower 5% and upper 5%
1117       if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,&paramT0,0.05,.95);
1118       fVMeanTimeSector[iSec]=paramT0[1];
1119     }
1120
1121     AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1122     AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1123     AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1124     AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1125
1126     TH2S *hQ   = GetHistoQ(iSec);
1127     TH2S *hRMS = GetHistoRMS(iSec);
1128
1129     Short_t *arrayhQ   = hQ->GetArray();
1130     Short_t *arrayhT0  = hT0->GetArray();
1131     Short_t *arrayhRMS = hRMS->GetArray();
1132
1133     UInt_t nChannels = fROC->GetNChannels(iSec);
1134     Float_t meanTsec = fVMeanTimeSector[iSec];
1135
1136   //debug
1137     Int_t row=0;
1138     Int_t pad=0;
1139     Int_t padc=0;
1140   //! debug
1141
1142     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1143
1144       Float_t cogTime0 = -1000;
1145       Float_t cogQ     = -1000;
1146       Float_t cogRMS   = -1000;
1147       Float_t cogOut   = 0;
1148
1149       Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1150       Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1151       Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1152 /*
1153       AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1154       AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1155       AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1156       cogQ     = paramQ[1];
1157       cogTime0 = paramT0[1];
1158       cogRMS   = paramRMS[1];
1159 */
1160       cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1161       cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1162       cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1163
1164       /*
1165       if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1166     cogOut = 1;
1167     cogTime0 = 0;
1168     cogQ     = 0;
1169     cogRMS   = 0;
1170       }
1171 */
1172 //       rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1173       rocQ->SetValue(iChannel, cogQ);
1174       rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1175       rocRMS->SetValue(iChannel, cogRMS);
1176       rocOut->SetValue(iChannel, cogOut);
1177       //debug
1178       if ( GetStreamLevel() > 2 ){
1179         TTreeSRedirector *streamer=GetDebugStreamer();
1180         if ( streamer ) {
1181           while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1182           pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1183           padc = pad-(fROC->GetNPads(iSec,row)/2);
1184
1185           (*streamer) << "DataEnd" <<
1186             "Sector="  << iSec      <<
1187             "Pad="     << pad       <<
1188             "PadC="    << padc      <<
1189             "Row="     << row       <<
1190             "PadSec="  << iChannel   <<
1191             "Q="       << cogQ      <<
1192             "T0="      << cogTime0  <<
1193             "RMS="     << cogRMS    <<
1194             "\n";
1195         }
1196       }
1197       //! debug
1198     }
1199
1200
1201   }
1202 }
1203 //_____________________________________________________________________
1204 //_________________________  Test Functions ___________________________
1205 //_____________________________________________________________________
1206 TObjArray* AliTPCCalibPulser::TestBinning()
1207 {
1208   //
1209   //  Function to test the binning of the reference histograms
1210   //  type: T0, Q or RMS
1211   //  mode: 0 - number of filled bins per channel
1212   //        1 - number of empty bins between filled bins in one ROC
1213   //  returns TObjArray with the test histograms type*2+mode:
1214   //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1215
1216
1217   TObjArray *histArray = new TObjArray(6);
1218   const Char_t *type[] = {"T0","Q","RMS"};
1219   Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1220
1221   for (Int_t itype = 0; itype<3; ++itype){
1222     for (Int_t imode=0; imode<2; ++imode){
1223       Int_t icount = itype*2+imode;
1224       histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1225                                 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1226                                 72,0,72),
1227                        icount);
1228     }
1229   }
1230
1231
1232   TH2S *hRef=0x0;
1233   Short_t *array=0x0;
1234   for (Int_t itype = 0; itype<3; ++itype){
1235     for (Int_t iSec=0; iSec<72; ++iSec){
1236       if ( itype == 0 ) hRef = GetHistoT0(iSec);
1237       if ( itype == 1 ) hRef = GetHistoQ(iSec);
1238       if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1239       if ( hRef == 0x0 ) continue;
1240       array = (hRef->GetArray());
1241       UInt_t nChannels = fROC->GetNChannels(iSec);
1242
1243       Int_t nempty=0;
1244       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1245         Int_t nfilled=0;
1246         Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1247         Int_t c1 = 0;
1248         Int_t c2 = 0;
1249         for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1250           if ( array[offset+iBin]>0 ) {
1251             nfilled++;
1252             if ( c1 && c2 ) nempty++;
1253             else c1 = 1;
1254           }
1255           else if ( c1 ) c2 = 1;
1256
1257         }
1258         ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1259       }
1260       ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1261     }
1262   }
1263   return histArray;
1264 }