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           "Sector=" <<fCurrentSector<<
 696           "Row="    <<fCurrentRow<<
 697           "Pad="    <<fCurrentPad<<
 698           "PadC="   <<padc<<
 699           "Channel="<<fCurrentChannel<<
 700           "Sum="    <<qSum<<
 701           "params.="<<&param<<
 702           "signal.=" <<&fPadSignal<<
 703           "\n";
 704       }
 705     } else { //debug > 1
 706       (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
 707       (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 708       (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
 709     }
 710   }
 711   ResetPad();
 712 }
 713 //_____________________________________________________________________
 714 void AliTPCCalibPulser::EndEvent()
 715 {
 716     //
 717     //  Process data of current event
 718     //
 719
 720     //check if last pad has allready been processed, if not do so
 721   if ( fMaxTimeBin>-1 ) ProcessPad();
 722
 723     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 724   for ( Int_t iSec = 0; iSec<72; ++iSec ){
 725     TVectorF *vTimes = GetPadTimesEvent(iSec);
 726     if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
 727     Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 728     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 729       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
 730
 731       GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
 732             //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
 733
 734
 735       //Debug start
 736       if ( GetStreamLevel()>1 ){
 737         TTreeSRedirector *streamer=GetDebugStreamer();
 738         if ( streamer ) {
 739           Int_t row=0;
 740           Int_t pad=0;
 741           Int_t padc=0;
 742
 743           Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 744           Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 745
 746           UInt_t channel=iChannel;
 747           Int_t sector=iSec;
 748
 749           while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 750           pad = channel-fROC->GetRowIndexes(sector)[row];
 751           padc = pad-(fROC->GetNPads(sector,row)/2);
 752
 753           (*streamer) << "DataPad" <<
 754 //                  "Event=" << fEvent <<
 755             "Sector="<< sector <<
 756             "Row="   << row<<
 757             "Pad="   << pad <<
 758             "PadC="  << padc <<
 759             "PadSec="<< channel <<
 760             "Time0="  << time0 <<
 761             "Time="  << time <<
 762             "RMS="   << rms <<
 763             "Sum="   << q <<
 764             "\n";
 765         }
 766       }
 767       //Debug end
 768     }
 769   }
 770 }
 771 //_____________________________________________________________________
 772 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
 773                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 774                                   const Char_t *type, Bool_t force)
 775 {
 776     //
 777     // return pointer to Q histogram
 778     // if force is true create a new histogram if it doesn't exist allready
 779     //
 780     if ( !force || arr->UncheckedAt(sector) )
 781         return (TH2S*)arr->UncheckedAt(sector);
 782
 783     // if we are forced and histogram doesn't yes exist create it
 784     Char_t name[255], title[255];
 785
 786     sprintf(name,"hCalib%s%.2d",type,sector);
 787     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 788
 789     // new histogram with Q calib information. One value for each pad!
 790     TH2S* hist = new TH2S(name,title,
 791                           nbinsY, ymin, ymax,
 792                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 793     hist->SetDirectory(0);
 794     arr->AddAt(hist,sector);
 795     return hist;
 796 }
 797 //_____________________________________________________________________
 798 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
 799 {
 800     //
 801     // return pointer to T0 histogram
 802     // if force is true create a new histogram if it doesn't exist allready
 803     //
 804     TObjArray *arr = &fHistoT0Array;
 805     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 806 }
 807 //_____________________________________________________________________
 808 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
 809 {
 810     //
 811     // return pointer to Q histogram
 812     // if force is true create a new histogram if it doesn't exist allready
 813     //
 814     TObjArray *arr = &fHistoQArray;
 815     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 816 }
 817 //_____________________________________________________________________
 818 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
 819 {
 820     //
 821     // return pointer to Q histogram
 822     // if force is true create a new histogram if it doesn't exist allready
 823     //
 824     TObjArray *arr = &fHistoRMSArray;
 825     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 826 }
 827 //_____________________________________________________________________
 828 TH2F* AliTPCCalibPulser::GetHistoTSec()
 829 {
 830     //
 831     // return the pointer to the abs time distribution per sector
 832     // create it if it does not exist
 833     //
 834     if ( !fHMeanTimeSector )   //!!!if you change the binning here, you should also change it in the Analyse Function!!
 835         fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
 836                                     20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
 837                                     72,0,72);
 838    return fHMeanTimeSector;
 839 }
 840 //_____________________________________________________________________
 841 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
 842 {
 843     //
 844     // return pointer to Pad Info from 'arr' for the current event and sector
 845     // if force is true create it if it doesn't exist allready
 846     //
 847     if ( !force || arr->UncheckedAt(sector) )
 848         return (TVectorF*)arr->UncheckedAt(sector);
 849
 850     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 851     arr->AddAt(vect,sector);
 852     return vect;
 853 }
 854 //_____________________________________________________________________
 855 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
 856 {
 857     //
 858     // return pointer to Pad Times Array for the current event and sector
 859     // if force is true create it if it doesn't exist allready
 860     //
 861     TObjArray *arr = &fPadTimesArrayEvent;
 862     return GetPadInfoEvent(sector,arr,force);
 863 }
 864 //_____________________________________________________________________
 865 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
 866 {
 867     //
 868     // return pointer to Pad Q Array for the current event and sector
 869     // if force is true create it if it doesn't exist allready
 870     // for debugging purposes only
 871     //
 872
 873     TObjArray *arr = &fPadQArrayEvent;
 874     return GetPadInfoEvent(sector,arr,force);
 875 }
 876 //_____________________________________________________________________
 877 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
 878 {
 879     //
 880     // return pointer to Pad RMS Array for the current event and sector
 881     // if force is true create it if it doesn't exist allready
 882     // for debugging purposes only
 883     //
 884     TObjArray *arr = &fPadRMSArrayEvent;
 885     return GetPadInfoEvent(sector,arr,force);
 886 }
 887 //_____________________________________________________________________
 888 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
 889 {
 890     //
 891     // return pointer to Pad RMS Array for the current event and sector
 892     // if force is true create it if it doesn't exist allready
 893     // for debugging purposes only
 894     //
 895     TObjArray *arr = &fPadPedestalArrayEvent;
 896     return GetPadInfoEvent(sector,arr,force);
 897 }
 898 //_____________________________________________________________________
 899 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
 900 {
 901     //
 902     // return pointer to ROC Calibration
 903     // if force is true create a new histogram if it doesn't exist allready
 904     //
 905     if ( !force || arr->UncheckedAt(sector) )
 906         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 907
 908     // if we are forced and histogram doesn't yes exist create it
 909
 910     // new AliTPCCalROC for T0 information. One value for each pad!
 911     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 912     arr->AddAt(croc,sector);
 913     return croc;
 914 }
 915 //_____________________________________________________________________
 916 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
 917 {
 918     //
 919     // return pointer to Carge ROC Calibration
 920     // if force is true create a new histogram if it doesn't exist allready
 921     //
 922     TObjArray *arr = &fCalRocArrayT0;
 923     return GetCalRoc(sector, arr, force);
 924 }
 925 //_____________________________________________________________________
 926 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
 927 {
 928     //
 929     // return pointer to T0 ROC Calibration
 930     // if force is true create a new histogram if it doesn't exist allready
 931     //
 932     TObjArray *arr = &fCalRocArrayQ;
 933     return GetCalRoc(sector, arr, force);
 934 }
 935 //_____________________________________________________________________
 936 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
 937 {
 938     //
 939     // return pointer to signal width ROC Calibration
 940     // if force is true create a new histogram if it doesn't exist allready
 941     //
 942     TObjArray *arr = &fCalRocArrayRMS;
 943     return GetCalRoc(sector, arr, force);
 944 }
 945 //_____________________________________________________________________
 946 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 947 {
 948     //
 949     // return pointer to Outliers
 950     // if force is true create a new histogram if it doesn't exist allready
 951     //
 952     TObjArray *arr = &fCalRocArrayOutliers;
 953     return GetCalRoc(sector, arr, force);
 954 }
 955 //_____________________________________________________________________
 956 void AliTPCCalibPulser::ResetEvent()
 957 {
 958     //
 959     //  Reset global counters  -- Should be called before each event is processed
 960     //
 961     fLastSector=-1;
 962     fCurrentSector=-1;
 963     fCurrentRow=-1;
 964     fCurrentPad=-1;
 965     fCurrentChannel=-1;
 966
 967     ResetPad();
 968
 969     fPadTimesArrayEvent.Delete();
 970
 971     fPadQArrayEvent.Delete();
 972     fPadRMSArrayEvent.Delete();
 973     fPadPedestalArrayEvent.Delete();
 974
 975     for ( Int_t i=0; i<72; ++i ){
 976         fVTime0Offset[i]=0;
 977         fVTime0OffsetCounter[i]=0;
 978     }
 979 }
 980 //_____________________________________________________________________
 981 void AliTPCCalibPulser::ResetPad()
 982 {
 983     //
 984     //  Reset pad infos -- Should be called after a pad has been processed
 985     //
 986     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 987         fPadSignal[i] = 0;
 988     fMaxTimeBin = -1;
 989     fMaxPadSignal = -1;
 990     fPadPedestal  = -1;
 991     fPadNoise     = -1;
 992 }
 993 //_____________________________________________________________________
 994 Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
 995 {
 996     //
 997     // return true if pad is on the edge of a row
 998     //
 999     Int_t edge1   = 0;
1000     Int_t edge2   = fROC->GetNPads(sector,row)-1;
1001     if ( pad == edge1 || pad == edge2 ) return kTRUE;
1002
1003     return kFALSE;
1004 }
1005 //_____________________________________________________________________
1006 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
1007 {
1008     //
1009     //  Merge reference histograms of sig to the current AliTPCCalibPulser
1010     //
1011
1012     //merge histograms
1013     for (Int_t iSec=0; iSec<72; ++iSec){
1014         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
1015         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
1016         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1017
1018
1019         if ( hRefQmerge ){
1020             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1021             TH2S *hRefQ   = GetHistoQ(iSec);
1022             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1023             else {
1024                 TH2S *hist = new TH2S(*hRefQmerge);
1025                 hist->SetDirectory(0);
1026                 fHistoQArray.AddAt(hist, iSec);
1027             }
1028             hRefQmerge->SetDirectory(dir);
1029         }
1030         if ( hRefT0merge ){
1031             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1032             TH2S *hRefT0  = GetHistoT0(iSec);
1033             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1034             else {
1035                 TH2S *hist = new TH2S(*hRefT0merge);
1036                 hist->SetDirectory(0);
1037                 fHistoT0Array.AddAt(hist, iSec);
1038             }
1039             hRefT0merge->SetDirectory(dir);
1040         }
1041         if ( hRefRMSmerge ){
1042             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1043             TH2S *hRefRMS = GetHistoRMS(iSec);
1044             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1045             else {
1046                 TH2S *hist = new TH2S(*hRefRMSmerge);
1047                 hist->SetDirectory(0);
1048                 fHistoRMSArray.AddAt(hist, iSec);
1049             }
1050             hRefRMSmerge->SetDirectory(dir);
1051         }
1052
1053     }
1054     if ( sig->fHMeanTimeSector ){
1055         TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1056         if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1057         else {
1058             fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1059             fHMeanTimeSector->SetDirectory(0);
1060         }
1061         sig->fHMeanTimeSector->SetDirectory(dir);
1062     }
1063 }
1064 //_____________________________________________________________________
1065 void AliTPCCalibPulser::Analyse()
1066 {
1067   //
1068   //  Calculate calibration constants
1069   //
1070
1071   TVectorD paramQ(3);
1072   TVectorD paramT0(3);
1073   TVectorD paramRMS(3);
1074   TMatrixD dummy(3,3);
1075   //calculate mean time for each sector and mean time for each side
1076   TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1077   fVMeanTimeSector.Zero();
1078
1079   for (Int_t iSec=0; iSec<72; ++iSec){
1080     TH2S *hT0 = GetHistoT0(iSec);
1081     if (!hT0 ) continue;
1082     //calculate sector mean T
1083     if ( fHMeanTimeSector ){
1084       Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1085       Int_t offset = (nbinsT+2)*(iSec+1);
1086       Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1087       Int_t entries=0;
1088       for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1089       hMeanTsec.Set(nbinsT+2,arrP);
1090       hMeanTsec.SetEntries(entries);
1091       paramT0.Zero();
1092       // truncated mean: remove lower 5% and upper 5%
1093       if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,&paramT0,0.05,.95);
1094       fVMeanTimeSector[iSec]=paramT0[1];
1095     }
1096
1097     AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1098     AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1099     AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1100     AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1101
1102     TH2S *hQ   = GetHistoQ(iSec);
1103     TH2S *hRMS = GetHistoRMS(iSec);
1104
1105     Short_t *arrayhQ   = hQ->GetArray();
1106     Short_t *arrayhT0  = hT0->GetArray();
1107     Short_t *arrayhRMS = hRMS->GetArray();
1108
1109     UInt_t nChannels = fROC->GetNChannels(iSec);
1110     Float_t meanTsec = fVMeanTimeSector[iSec];
1111
1112   //debug
1113     Int_t row=0;
1114     Int_t pad=0;
1115     Int_t padc=0;
1116   //! debug
1117
1118     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1119
1120       Float_t cogTime0 = -1000;
1121       Float_t cogQ     = -1000;
1122       Float_t cogRMS   = -1000;
1123       Float_t cogOut   = 0;
1124
1125       Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1126       Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1127       Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1128 /*
1129       AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1130       AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1131       AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1132       cogQ     = paramQ[1];
1133       cogTime0 = paramT0[1];
1134       cogRMS   = paramRMS[1];
1135 */
1136       cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1137       cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1138       cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1139
1140       /*
1141       if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1142     cogOut = 1;
1143     cogTime0 = 0;
1144     cogQ     = 0;
1145     cogRMS   = 0;
1146       }
1147 */
1148 //       rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1149       rocQ->SetValue(iChannel, cogQ);
1150       rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1151       rocRMS->SetValue(iChannel, cogRMS);
1152       rocOut->SetValue(iChannel, cogOut);
1153       //debug
1154       if ( GetStreamLevel() > 2 ){
1155         TTreeSRedirector *streamer=GetDebugStreamer();
1156         if ( streamer ) {
1157           while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1158           pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1159           padc = pad-(fROC->GetNPads(iSec,row)/2);
1160
1161           (*streamer) << "DataEnd" <<
1162             "Sector="  << iSec      <<
1163             "Pad="     << pad       <<
1164             "PadC="    << padc      <<
1165             "Row="     << row       <<
1166             "PadSec="  << iChannel   <<
1167             "Q="       << cogQ      <<
1168             "T0="      << cogTime0  <<
1169             "RMS="     << cogRMS    <<
1170             "\n";
1171         }
1172       }
1173       //! debug
1174     }
1175
1176
1177   }
1178 }
1179 //_____________________________________________________________________
1180 //_________________________  Test Functions ___________________________
1181 //_____________________________________________________________________
1182 TObjArray* AliTPCCalibPulser::TestBinning()
1183 {
1184   //
1185   //  Function to test the binning of the reference histograms
1186   //  type: T0, Q or RMS
1187   //  mode: 0 - number of filled bins per channel
1188   //        1 - number of empty bins between filled bins in one ROC
1189   //  returns TObjArray with the test histograms type*2+mode:
1190   //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1191
1192
1193   TObjArray *histArray = new TObjArray(6);
1194   const Char_t *type[] = {"T0","Q","RMS"};
1195   Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1196
1197   for (Int_t itype = 0; itype<3; ++itype){
1198     for (Int_t imode=0; imode<2; ++imode){
1199       Int_t icount = itype*2+imode;
1200       histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1201                                 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1202                                 72,0,72),
1203                        icount);
1204     }
1205   }
1206
1207
1208   TH2S *hRef=0x0;
1209   Short_t *array=0x0;
1210   for (Int_t itype = 0; itype<3; ++itype){
1211     for (Int_t iSec=0; iSec<72; ++iSec){
1212       if ( itype == 0 ) hRef = GetHistoT0(iSec);
1213       if ( itype == 1 ) hRef = GetHistoQ(iSec);
1214       if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1215       if ( hRef == 0x0 ) continue;
1216       array = (hRef->GetArray());
1217       UInt_t nChannels = fROC->GetNChannels(iSec);
1218
1219       Int_t nempty=0;
1220       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1221         Int_t nfilled=0;
1222         Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1223         Int_t c1 = 0;
1224         Int_t c2 = 0;
1225         for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1226           if ( array[offset+iBin]>0 ) {
1227             nfilled++;
1228             if ( c1 && c2 ) nempty++;
1229             else c1 = 1;
1230           }
1231           else if ( c1 ) c2 = 1;
1232
1233         }
1234         ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1235       }
1236       ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1237     }
1238   }
1239   return histArray;
1240 }