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 "AliTPCCalROC.h"
 186 #include "AliTPCCalPad.h"
 187 #include "AliTPCROC.h"
 188 #include "AliTPCParam.h"
 189 #include "AliTPCCalibPulser.h"
 190 #include "AliTPCcalibDB.h"
 191 #include "AliMathBase.h"
 192 #include "AliLog.h"
 193 #include "TTreeStream.h"
 194
 195 //date
 196 #include "event.h"
 197
 198
 199
 200
 201 ClassImp(AliTPCCalibPulser)
 202
 203 AliTPCCalibPulser::AliTPCCalibPulser() :
 204 AliTPCCalibRawBase(),
 205 fNbinsT0(200),
 206 fXminT0(-2),
 207 fXmaxT0(2),
 208 fNbinsQ(200),
 209 fXminQ(10),
 210 fXmaxQ(40),
 211 fNbinsRMS(100),
 212 fXminRMS(0.1),
 213 fXmaxRMS(5.1),
 214 fPeakIntMinus(2),
 215 fPeakIntPlus(2),
 216 fIsZeroSuppressed(kFALSE),
 217 fLastSector(-1),
 218 fParam(new AliTPCParam),
 219 fPedestalTPC(0x0),
 220 fPadNoiseTPC(0x0),
 221 fOutliers(0x0),
 222 fPedestalROC(0x0),
 223 fPadNoiseROC(0x0),
 224 fCalRocArrayT0(72),
 225 fCalRocArrayQ(72),
 226 fCalRocArrayRMS(72),
 227 fCalRocArrayOutliers(72),
 228 fHistoQArray(72),
 229 fHistoT0Array(72),
 230 fHistoRMSArray(72),
 231 fHMeanTimeSector(0x0),
 232 fVMeanTimeSector(72),
 233 fPadTimesArrayEvent(72),
 234 fPadQArrayEvent(72),
 235 fPadRMSArrayEvent(72),
 236 fPadPedestalArrayEvent(72),
 237 fCurrentChannel(-1),
 238 fCurrentSector(-1),
 239 fCurrentRow(-1),
 240 fCurrentPad(-1),
 241 fMaxPadSignal(-1),
 242 fMaxTimeBin(-1),
 243 fPadSignal(1024),
 244 fPadPedestal(0),
 245 fPadNoise(0),
 246 fVTime0Offset(72),
 247 fVTime0OffsetCounter(72)
 248 {
 249   //
 250   // AliTPCSignal default constructor
 251   //
 252   SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
 253   fFirstTimeBin=60;
 254   fLastTimeBin=900;
 255   fParam->Update();
 256 }
 257 //_____________________________________________________________________
 258 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
 259 AliTPCCalibRawBase(sig),
 260 fNbinsT0(sig.fNbinsT0),
 261 fXminT0(sig.fXminT0),
 262 fXmaxT0(sig.fXmaxT0),
 263 fNbinsQ(sig.fNbinsQ),
 264 fXminQ(sig.fXminQ),
 265 fXmaxQ(sig.fXmaxQ),
 266 fNbinsRMS(sig.fNbinsRMS),
 267 fXminRMS(sig.fXminRMS),
 268 fXmaxRMS(sig.fXmaxRMS),
 269 fPeakIntMinus(sig.fPeakIntMinus),
 270 fPeakIntPlus(sig.fPeakIntPlus),
 271 fIsZeroSuppressed(sig.fIsZeroSuppressed),
 272 fLastSector(-1),
 273 fParam(new AliTPCParam),
 274 fPedestalTPC(0x0),
 275 fPadNoiseTPC(0x0),
 276 fOutliers(0x0),
 277 fPedestalROC(0x0),
 278 fPadNoiseROC(0x0),
 279 fCalRocArrayT0(72),
 280 fCalRocArrayQ(72),
 281 fCalRocArrayRMS(72),
 282 fCalRocArrayOutliers(72),
 283 fHistoQArray(72),
 284 fHistoT0Array(72),
 285 fHistoRMSArray(72),
 286 fHMeanTimeSector(0x0),
 287 fVMeanTimeSector(72),
 288 fPadTimesArrayEvent(72),
 289 fPadQArrayEvent(72),
 290 fPadRMSArrayEvent(72),
 291 fPadPedestalArrayEvent(72),
 292 fCurrentChannel(-1),
 293 fCurrentSector(-1),
 294 fCurrentRow(-1),
 295 fCurrentPad(-1),
 296 fMaxPadSignal(-1),
 297 fMaxTimeBin(-1),
 298 fPadSignal(1024),
 299 fPadPedestal(0),
 300 fPadNoise(0),
 301 fVTime0Offset(72),
 302 fVTime0OffsetCounter(72)
 303 {
 304   //
 305   // AliTPCSignal default constructor
 306   //
 307
 308   for (Int_t iSec = 0; iSec < 72; ++iSec){
 309     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 310     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 311     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 312     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 313
 314     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 315     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 316     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 317
 318     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 319     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 320     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 321     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 322
 323     if ( hQ != 0x0 ){
 324       TH2S *hNew = new TH2S(*hQ);
 325       hNew->SetDirectory(0);
 326       fHistoQArray.AddAt(hNew,iSec);
 327     }
 328     if ( hT0 != 0x0 ){
 329       TH2S *hNew = new TH2S(*hT0);
 330       hNew->SetDirectory(0);
 331       fHistoQArray.AddAt(hNew,iSec);
 332     }
 333     if ( hRMS != 0x0 ){
 334       TH2S *hNew = new TH2S(*hRMS);
 335       hNew->SetDirectory(0);
 336       fHistoQArray.AddAt(hNew,iSec);
 337     }
 338     fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
 339   }
 340
 341   if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
 342   fParam->Update();
 343 }
 344 AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
 345 AliTPCCalibRawBase(),
 346 fNbinsT0(200),
 347 fXminT0(-2),
 348 fXmaxT0(2),
 349 fNbinsQ(200),
 350 fXminQ(10),
 351 fXmaxQ(40),
 352 fNbinsRMS(100),
 353 fXminRMS(0.1),
 354 fXmaxRMS(5.1),
 355 fPeakIntMinus(2),
 356 fPeakIntPlus(2),
 357 fIsZeroSuppressed(kFALSE),
 358 fLastSector(-1),
 359 fParam(new  AliTPCParam),
 360 fPedestalTPC(0x0),
 361 fPadNoiseTPC(0x0),
 362 fOutliers(0x0),
 363 fPedestalROC(0x0),
 364 fPadNoiseROC(0x0),
 365 fCalRocArrayT0(72),
 366 fCalRocArrayQ(72),
 367 fCalRocArrayRMS(72),
 368 fCalRocArrayOutliers(72),
 369 fHistoQArray(72),
 370 fHistoT0Array(72),
 371 fHistoRMSArray(72),
 372 fHMeanTimeSector(0x0),
 373 fVMeanTimeSector(72),
 374 fPadTimesArrayEvent(72),
 375 fPadQArrayEvent(72),
 376 fPadRMSArrayEvent(72),
 377 fPadPedestalArrayEvent(72),
 378 fCurrentChannel(-1),
 379 fCurrentSector(-1),
 380 fCurrentRow(-1),
 381 fCurrentPad(-1),
 382 fMaxPadSignal(-1),
 383 fMaxTimeBin(-1),
 384 fPadSignal(1024),
 385 fPadPedestal(0),
 386 fPadNoise(0),
 387 fVTime0Offset(72),
 388 fVTime0OffsetCounter(72)
 389 {
 390   //
 391   // This constructor uses a TMap for setting some parametes
 392   //
 393   SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
 394   fFirstTimeBin=60;
 395   fLastTimeBin=900;
 396   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 397   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 398   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 399   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 400   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 401   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 402   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 403   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 404   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 405   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 406   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 407   if (config->GetValue("PeakIntMinus")) fPeakIntMinus = (Int_t)((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atof();
 408   if (config->GetValue("PeakIntPlus")) fPeakIntPlus = (Int_t)((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atof();
 409   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 410
 411   fParam->Update();
 412 }
 413 //_____________________________________________________________________
 414 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
 415 {
 416   //
 417   // assignment operator
 418   //
 419   if (&source == this) return *this;
 420   new (this) AliTPCCalibPulser(source);
 421
 422   return *this;
 423 }
 424 //_____________________________________________________________________
 425 AliTPCCalibPulser::~AliTPCCalibPulser()
 426 {
 427   //
 428   // destructor
 429   //
 430
 431   Reset();
 432   delete fParam;
 433 }
 434 void AliTPCCalibPulser::Reset()
 435 {
 436   //
 437   // Delete all information: Arrays, Histograms, CalRoc objects
 438   //
 439   fCalRocArrayT0.Delete();
 440   fCalRocArrayQ.Delete();
 441   fCalRocArrayRMS.Delete();
 442   fCalRocArrayOutliers.Delete();
 443
 444   fHistoQArray.Delete();
 445   fHistoT0Array.Delete();
 446   fHistoRMSArray.Delete();
 447
 448   fPadTimesArrayEvent.Delete();
 449   fPadQArrayEvent.Delete();
 450   fPadRMSArrayEvent.Delete();
 451   fPadPedestalArrayEvent.Delete();
 452
 453   if (fHMeanTimeSector) delete fHMeanTimeSector;
 454 }
 455 //_____________________________________________________________________
 456 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
 457                                 const Int_t icRow,
 458                                 const Int_t icPad,
 459                                 const Int_t icTimeBin,
 460                                 const Float_t csignal)
 461 {
 462     //
 463     // Signal filling methode on the fly pedestal and time offset correction if necessary.
 464     // no extra analysis necessary. Assumes knowledge of the signal shape!
 465     // assumes that it is looped over consecutive time bins of one pad
 466     //
 467
 468   if (icRow<0) return 0;
 469   if (icPad<0) return 0;
 470   if (icTimeBin<0) return 0;
 471   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 472
 473   if ( icRow<0 || icPad<0 ){
 474     AliWarning("Wrong Pad or Row number, skipping!");
 475     return 0;
 476   }
 477
 478   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 479
 480     //init first pad and sector in this event
 481   if ( fCurrentChannel == -1 ) {
 482     fCurrentChannel = iChannel;
 483     fCurrentSector  = icsector;
 484     fCurrentRow     = icRow;
 485     fCurrentPad     = icPad;
 486   }
 487
 488     //process last pad if we change to a new one
 489   if ( iChannel != fCurrentChannel ){
 490     ProcessPad();
 491     fLastSector=fCurrentSector;
 492     fCurrentChannel = iChannel;
 493     fCurrentSector  = icsector;
 494     fCurrentRow     = icRow;
 495     fCurrentPad     = icPad;
 496   }
 497
 498     //fill signals for current pad
 499   fPadSignal[icTimeBin]=csignal;
 500   if ( csignal > fMaxPadSignal ){
 501     fMaxPadSignal = csignal;
 502     fMaxTimeBin   = icTimeBin;
 503   }
 504   return 0;
 505 }
 506 //_____________________________________________________________________
 507 void AliTPCCalibPulser::FindPedestal(Float_t part)
 508 {
 509     //
 510     // find pedestal and noise for the current pad. Use either database or
 511     // truncated mean with part*100%
 512     //
 513   Bool_t noPedestal = kTRUE;;
 514   if (fPedestalTPC&&fPadNoiseTPC){
 515         //use pedestal database
 516         //only load new pedestals if the sector has changed
 517     if ( fCurrentSector!=fLastSector ){
 518       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 519       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 520     }
 521
 522     if ( fPedestalROC&&fPadNoiseROC ){
 523       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 524       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 525       noPedestal   = kFALSE;
 526     }
 527
 528   }
 529
 530     //if we are not running with pedestal database, or for the current sector there is no information
 531     //available, calculate the pedestal and noise on the fly
 532   if ( noPedestal ) {
 533     const Int_t kPedMax = 100;  //maximum pedestal value
 534     Float_t  max    =  0;
 535     Float_t  maxPos =  0;
 536     Int_t    median =  -1;
 537     Int_t    count0 =  0;
 538     Int_t    count1 =  0;
 539   //
 540     Float_t padSignal=0;
 541         //
 542     UShort_t histo[kPedMax];
 543     memset(histo,0,kPedMax*sizeof(UShort_t));
 544
 545     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 546       padSignal = fPadSignal.GetMatrixArray()[i];
 547       if (padSignal<=0) continue;
 548       if (padSignal>max && i>10) {
 549         max = padSignal;
 550         maxPos = i;
 551       }
 552       if (padSignal>kPedMax-1) continue;
 553       histo[Int_t(padSignal+0.5)]++;
 554       count0++;
 555     }
 556       //
 557     for (Int_t i=1; i<kPedMax; ++i){
 558       if (count1<count0*0.5) median=i;
 559       count1+=histo[i];
 560     }
 561   // truncated mean
 562   //
 563         // what if by chance histo[median] == 0 ?!?
 564     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 565   //
 566     for (Int_t idelta=1; idelta<10; ++idelta){
 567       if (median-idelta<=0) continue;
 568       if (median+idelta>kPedMax) continue;
 569       if (count<part*count1){
 570         count+=histo[median-idelta];
 571         mean +=histo[median-idelta]*(median-idelta);
 572         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 573         count+=histo[median+idelta];
 574         mean +=histo[median+idelta]*(median+idelta);
 575         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 576       }
 577     }
 578     fPadPedestal = 0;
 579     fPadNoise    = 0;
 580     if ( count > 0 ) {
 581       mean/=count;
 582       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 583       fPadPedestal = mean;
 584       fPadNoise    = rms;
 585     }
 586   }
 587   fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
 588 }
 589 //_____________________________________________________________________
 590 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
 591 {
 592 //
 593     //  Find position, signal width and height of the CE signal (last signal)
 594     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 595     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 596     //
 597
 598   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 599   Int_t   cemaxpos       = fMaxTimeBin;
 600   Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal sum
 601   Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal max
 602   const Int_t    kCemin  = fPeakIntMinus;             // range for the analysis of the ce signal +- channels from the peak
 603   const Int_t    kCemax  = fPeakIntPlus;
 604   param[0] = ceQmax;
 605   param[1] = ceTime;
 606   param[2] = ceRMS;
 607   qSum     = ceQsum;
 608
 609   if (cemaxpos>0){
 610     ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 611     if ( ceQmax<ceMaxThreshold ) return;
 612     for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
 613       Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 614       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 615         ceTime+=signal*(i+0.5);
 616         ceRMS +=signal*(i+0.5)*(i+0.5);
 617         ceQsum+=signal;
 618       }
 619     }
 620   }
 621   if (ceQsum>ceSumThreshold) {
 622     ceTime/=ceQsum;
 623     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 624     ceTime-=GetL1PhaseTB();
 625         //only fill the Time0Offset if pad was not marked as an outlier!
 626     if ( !fOutliers ){
 627       //skip edge pads for calculating the mean time
 628       if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
 629         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 630         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 631         GetHistoTSec()->Fill(ceTime,fCurrentSector);
 632       }
 633     } else {
 634       if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
 635         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 636         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 637       }
 638     }
 639
 640   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
 641   //                                the pick-up signal should scale with the pad area. In addition
 642   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
 643   //                                ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
 644   //                                to the OROC signal (factor 2/3 for the IROCs).
 645     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
 646     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
 647
 648     ceQsum/=norm;
 649   } else {
 650     ceQmax=0;
 651     ceTime=0;
 652     ceRMS =0;
 653     ceQsum=0;
 654   }
 655   param[0] = ceQmax;
 656   param[1] = ceTime;
 657   param[2] = ceRMS;
 658   qSum     = ceQsum;
 659 }
 660 //_____________________________________________________________________
 661 void AliTPCCalibPulser::ProcessPad()
 662 {
 663   //
 664   //  Process data of current pad
 665   //
 666
 667   FindPedestal();
 668   TVectorD param(3);
 669   Float_t  qSum;
 670   FindPulserSignal(param, qSum);
 671
 672   Double_t meanT  = param[1];
 673   Double_t sigmaT = param[2];
 674
 675
 676   //Fill Event T0 counter
 677   (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 678
 679   //Fill Q histogram
 680 //   GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
 681   GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
 682
 683   //Fill RMS histogram
 684   GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 685
 686
 687     //Fill debugging info
 688   if ( GetStreamLevel()>0 ){
 689     TTreeSRedirector *streamer=GetDebugStreamer();
 690     if ( GetStreamLevel() == 1 ){
 691       if ( streamer ) {
 692         Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
 693         (*streamer) << "PadSignals" <<
 694           "Event="  <<fNevents <<
 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=" << fNevents <<
 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   ++fNevents;
 771 }
 772 //_____________________________________________________________________
 773 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
 774                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 775                                   const Char_t *type, Bool_t force)
 776 {
 777     //
 778     // return pointer to Q histogram
 779     // if force is true create a new histogram if it doesn't exist allready
 780     //
 781   if ( !force || arr->UncheckedAt(sector) )
 782     return (TH2S*)arr->UncheckedAt(sector);
 783
 784   // if we are forced and histogram doesn't yes exist create it
 785   // new histogram with Q calib information. One value for each pad!
 786   TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
 787                         nbinsY, ymin, ymax,
 788                         fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 789   hist->SetDirectory(0);
 790   arr->AddAt(hist,sector);
 791   return hist;
 792 }
 793 //_____________________________________________________________________
 794 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
 795 {
 796     //
 797     // return pointer to T0 histogram
 798     // if force is true create a new histogram if it doesn't exist allready
 799     //
 800   TObjArray *arr = &fHistoT0Array;
 801   return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 802 }
 803 //_____________________________________________________________________
 804 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
 805 {
 806     //
 807     // return pointer to Q histogram
 808     // if force is true create a new histogram if it doesn't exist allready
 809     //
 810   TObjArray *arr = &fHistoQArray;
 811   return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 812 }
 813 //_____________________________________________________________________
 814 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
 815 {
 816     //
 817     // return pointer to Q histogram
 818     // if force is true create a new histogram if it doesn't exist allready
 819     //
 820   TObjArray *arr = &fHistoRMSArray;
 821   return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 822 }
 823 //_____________________________________________________________________
 824 TH2F* AliTPCCalibPulser::GetHistoTSec()
 825 {
 826     //
 827     // return the pointer to the abs time distribution per sector
 828     // create it if it does not exist
 829     //
 830   if ( !fHMeanTimeSector )   //!!!if you change the binning here, you should also change it in the Analyse Function!!
 831     fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
 832                                 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
 833                                 72,0,72);
 834   return fHMeanTimeSector;
 835 }
 836 //_____________________________________________________________________
 837 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
 838 {
 839     //
 840     // return pointer to Pad Info from 'arr' for the current event and sector
 841     // if force is true create it if it doesn't exist allready
 842     //
 843   if ( !force || arr->UncheckedAt(sector) )
 844     return (TVectorF*)arr->UncheckedAt(sector);
 845
 846   TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 847   arr->AddAt(vect,sector);
 848   return vect;
 849 }
 850 //_____________________________________________________________________
 851 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
 852 {
 853     //
 854     // return pointer to Pad Times Array for the current event and sector
 855     // if force is true create it if it doesn't exist allready
 856     //
 857   TObjArray *arr = &fPadTimesArrayEvent;
 858   return GetPadInfoEvent(sector,arr,force);
 859 }
 860 //_____________________________________________________________________
 861 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
 862 {
 863     //
 864     // return pointer to Pad Q Array for the current event and sector
 865     // if force is true create it if it doesn't exist allready
 866     // for debugging purposes only
 867     //
 868
 869   TObjArray *arr = &fPadQArrayEvent;
 870   return GetPadInfoEvent(sector,arr,force);
 871 }
 872 //_____________________________________________________________________
 873 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
 874 {
 875     //
 876     // return pointer to Pad RMS Array for the current event and sector
 877     // if force is true create it if it doesn't exist allready
 878     // for debugging purposes only
 879     //
 880   TObjArray *arr = &fPadRMSArrayEvent;
 881   return GetPadInfoEvent(sector,arr,force);
 882 }
 883 //_____________________________________________________________________
 884 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
 885 {
 886     //
 887     // return pointer to Pad RMS Array for the current event and sector
 888     // if force is true create it if it doesn't exist allready
 889     // for debugging purposes only
 890     //
 891   TObjArray *arr = &fPadPedestalArrayEvent;
 892   return GetPadInfoEvent(sector,arr,force);
 893 }
 894 //_____________________________________________________________________
 895 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
 896 {
 897     //
 898     // return pointer to ROC Calibration
 899     // if force is true create a new histogram if it doesn't exist allready
 900     //
 901   if ( !force || arr->UncheckedAt(sector) )
 902     return (AliTPCCalROC*)arr->UncheckedAt(sector);
 903
 904     // if we are forced and histogram doesn't yes exist create it
 905
 906     // new AliTPCCalROC for T0 information. One value for each pad!
 907   AliTPCCalROC *croc = new AliTPCCalROC(sector);
 908   arr->AddAt(croc,sector);
 909   return croc;
 910 }
 911 //_____________________________________________________________________
 912 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
 913 {
 914     //
 915     // return pointer to Carge ROC Calibration
 916     // if force is true create a new histogram if it doesn't exist allready
 917     //
 918   TObjArray *arr = &fCalRocArrayT0;
 919   return GetCalRoc(sector, arr, force);
 920 }
 921 //_____________________________________________________________________
 922 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
 923 {
 924     //
 925     // return pointer to T0 ROC Calibration
 926     // if force is true create a new histogram if it doesn't exist allready
 927     //
 928   TObjArray *arr = &fCalRocArrayQ;
 929   return GetCalRoc(sector, arr, force);
 930 }
 931 //_____________________________________________________________________
 932 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
 933 {
 934     //
 935     // return pointer to signal width ROC Calibration
 936     // if force is true create a new histogram if it doesn't exist allready
 937     //
 938   TObjArray *arr = &fCalRocArrayRMS;
 939   return GetCalRoc(sector, arr, force);
 940 }
 941 //_____________________________________________________________________
 942 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 943 {
 944     //
 945     // return pointer to Outliers
 946     // if force is true create a new histogram if it doesn't exist allready
 947     //
 948   TObjArray *arr = &fCalRocArrayOutliers;
 949   return GetCalRoc(sector, arr, force);
 950 }
 951 //_____________________________________________________________________
 952 void AliTPCCalibPulser::ResetEvent()
 953 {
 954     //
 955     //  Reset global counters  -- Should be called before each event is processed
 956     //
 957   fLastSector=-1;
 958   fCurrentSector=-1;
 959   fCurrentRow=-1;
 960   fCurrentPad=-1;
 961   fCurrentChannel=-1;
 962
 963   ResetPad();
 964
 965   fPadTimesArrayEvent.Delete();
 966
 967   fPadQArrayEvent.Delete();
 968   fPadRMSArrayEvent.Delete();
 969   fPadPedestalArrayEvent.Delete();
 970
 971   for ( Int_t i=0; i<72; ++i ){
 972     fVTime0Offset[i]=0;
 973     fVTime0OffsetCounter[i]=0;
 974   }
 975 }
 976 //_____________________________________________________________________
 977 void AliTPCCalibPulser::ResetPad()
 978 {
 979     //
 980     //  Reset pad infos -- Should be called after a pad has been processed
 981     //
 982   for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 983     fPadSignal[i] = 0;
 984   fMaxTimeBin = -1;
 985   fMaxPadSignal = -1;
 986   fPadPedestal  = -1;
 987   fPadNoise     = -1;
 988 }
 989 //_____________________________________________________________________
 990 Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
 991 {
 992     //
 993     // return true if pad is on the edge of a row
 994     //
 995   Int_t edge1   = 0;
 996   Int_t edge2   = fROC->GetNPads(sector,row)-1;
 997   if ( pad == edge1 || pad == edge2 ) return kTRUE;
 998
 999   return kFALSE;
1000 }
1001 //_____________________________________________________________________
1002 void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
1003 {
1004   //
1005   //  Merge reference histograms of sig to the current AliTPCCalibPulser
1006   //
1007
1008   MergeBase(sig);
1009   //merge histograms
1010   for (Int_t iSec=0; iSec<72; ++iSec){
1011     TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
1012     TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
1013     TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1014
1015
1016     if ( hRefQmerge ){
1017       TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1018       TH2S *hRefQ   = GetHistoQ(iSec);
1019       if ( hRefQ ) hRefQ->Add(hRefQmerge);
1020       else {
1021         TH2S *hist = new TH2S(*hRefQmerge);
1022         hist->SetDirectory(0);
1023         fHistoQArray.AddAt(hist, iSec);
1024       }
1025       hRefQmerge->SetDirectory(dir);
1026     }
1027     if ( hRefT0merge ){
1028       TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1029       TH2S *hRefT0  = GetHistoT0(iSec);
1030       if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1031       else {
1032         TH2S *hist = new TH2S(*hRefT0merge);
1033         hist->SetDirectory(0);
1034         fHistoT0Array.AddAt(hist, iSec);
1035       }
1036       hRefT0merge->SetDirectory(dir);
1037     }
1038     if ( hRefRMSmerge ){
1039       TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1040       TH2S *hRefRMS = GetHistoRMS(iSec);
1041       if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1042       else {
1043         TH2S *hist = new TH2S(*hRefRMSmerge);
1044         hist->SetDirectory(0);
1045         fHistoRMSArray.AddAt(hist, iSec);
1046       }
1047       hRefRMSmerge->SetDirectory(dir);
1048     }
1049
1050   }
1051   if ( sig->fHMeanTimeSector ){
1052     TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1053     if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1054     else {
1055       fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1056       fHMeanTimeSector->SetDirectory(0);
1057     }
1058     sig->fHMeanTimeSector->SetDirectory(dir);
1059   }
1060 }
1061
1062
1063 //_____________________________________________________________________
1064 Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1065 {
1066   //
1067   // Merge all objects of this type in list
1068   //
1069
1070   Long64_t nmerged=1;
1071
1072   TIter next(list);
1073   AliTPCCalibPulser *ce=0;
1074   TObject *o=0;
1075
1076   while ( (o=next()) ){
1077     ce=dynamic_cast<AliTPCCalibPulser*>(o);
1078     if (ce){
1079       Merge(ce);
1080       ++nmerged;
1081     }
1082   }
1083
1084   return nmerged;
1085 }
1086
1087 //_____________________________________________________________________
1088 void AliTPCCalibPulser::Analyse()
1089 {
1090   //
1091   //  Calculate calibration constants
1092   //
1093
1094   TVectorD paramQ(3);
1095   TVectorD paramT0(3);
1096   TVectorD paramRMS(3);
1097   TMatrixD dummy(3,3);
1098   //calculate mean time for each sector and mean time for each side
1099   TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1100   fVMeanTimeSector.Zero();
1101
1102   for (Int_t iSec=0; iSec<72; ++iSec){
1103     TH2S *hT0 = GetHistoT0(iSec);
1104     if (!hT0 ) continue;
1105     //calculate sector mean T
1106     if ( fHMeanTimeSector ){
1107       Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1108       Int_t offset = (nbinsT+2)*(iSec+1);
1109       Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1110       Int_t entries=0;
1111       for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1112       hMeanTsec.Set(nbinsT+2,arrP);
1113       hMeanTsec.SetEntries(entries);
1114       paramT0.Zero();
1115       // truncated mean: remove lower 5% and upper 5%
1116       if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,&paramT0,0.05,.95);
1117       fVMeanTimeSector[iSec]=paramT0[1];
1118     }
1119
1120     AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1121     AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1122     AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1123     AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1124
1125     TH2S *hQ   = GetHistoQ(iSec);
1126     TH2S *hRMS = GetHistoRMS(iSec);
1127
1128     Short_t *arrayhQ   = hQ->GetArray();
1129     Short_t *arrayhT0  = hT0->GetArray();
1130     Short_t *arrayhRMS = hRMS->GetArray();
1131
1132     UInt_t nChannels = fROC->GetNChannels(iSec);
1133     Float_t meanTsec = fVMeanTimeSector[iSec];
1134
1135   //debug
1136     Int_t row=0;
1137     Int_t pad=0;
1138     Int_t padc=0;
1139   //! debug
1140
1141     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1142
1143       Float_t cogTime0 = -1000;
1144       Float_t cogQ     = -1000;
1145       Float_t cogRMS   = -1000;
1146       Float_t cogOut   = 0;
1147
1148       Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1149       Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1150       Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1151 /*
1152       AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1153       AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1154       AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1155       cogQ     = paramQ[1];
1156       cogTime0 = paramT0[1];
1157       cogRMS   = paramRMS[1];
1158 */
1159       cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1160       cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1161       cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1162
1163       /*
1164       if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1165     cogOut = 1;
1166     cogTime0 = 0;
1167     cogQ     = 0;
1168     cogRMS   = 0;
1169       }
1170 */
1171 //       rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1172       rocQ->SetValue(iChannel, cogQ);
1173       rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1174       rocRMS->SetValue(iChannel, cogRMS);
1175       rocOut->SetValue(iChannel, cogOut);
1176
1177       // in case a channel has no data set the value to 0
1178       if (TMath::Abs(cogTime0-fXminT0)<1e-10){
1179         rocQ->SetValue(iChannel, 0);
1180         rocT0->SetValue(iChannel, 0); //offset by mean time of the sector
1181         rocRMS->SetValue(iChannel, 0);
1182       }
1183
1184       //debug
1185       if ( GetStreamLevel() > 2 ){
1186         TTreeSRedirector *streamer=GetDebugStreamer();
1187         if ( streamer ) {
1188           while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1189           pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1190           padc = pad-(fROC->GetNPads(iSec,row)/2);
1191
1192           (*streamer) << "DataEnd" <<
1193             "Sector="  << iSec      <<
1194             "Pad="     << pad       <<
1195             "PadC="    << padc      <<
1196             "Row="     << row       <<
1197             "PadSec="  << iChannel   <<
1198             "Q="       << cogQ      <<
1199             "T0="      << cogTime0  <<
1200             "RMS="     << cogRMS    <<
1201             "\n";
1202         }
1203       }
1204       //! debug
1205     }
1206
1207
1208   }
1209 }
1210 //_____________________________________________________________________
1211 //_________________________  Test Functions ___________________________
1212 //_____________________________________________________________________
1213 TObjArray* AliTPCCalibPulser::TestBinning()
1214 {
1215   //
1216   //  Function to test the binning of the reference histograms
1217   //  type: T0, Q or RMS
1218   //  mode: 0 - number of filled bins per channel
1219   //        1 - number of empty bins between filled bins in one ROC
1220   //  returns TObjArray with the test histograms type*2+mode:
1221   //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1222
1223
1224   TObjArray *histArray = new TObjArray(6);
1225   const Char_t *type[] = {"T0","Q","RMS"};
1226   Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1227
1228   for (Int_t itype = 0; itype<3; ++itype){
1229     for (Int_t imode=0; imode<2; ++imode){
1230       Int_t icount = itype*2+imode;
1231       histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1232                                 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1233                                 72,0,72),
1234                        icount);
1235     }
1236   }
1237
1238
1239   TH2S *hRef=0x0;
1240   Short_t *array=0x0;
1241   for (Int_t itype = 0; itype<3; ++itype){
1242     for (Int_t iSec=0; iSec<72; ++iSec){
1243       if ( itype == 0 ) hRef = GetHistoT0(iSec);
1244       if ( itype == 1 ) hRef = GetHistoQ(iSec);
1245       if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1246       if ( hRef == 0x0 ) continue;
1247       array = (hRef->GetArray());
1248       UInt_t nChannels = fROC->GetNChannels(iSec);
1249
1250       Int_t nempty=0;
1251       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1252         Int_t nfilled=0;
1253         Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1254         Int_t c1 = 0;
1255         Int_t c2 = 0;
1256         for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1257           if ( array[offset+iBin]>0 ) {
1258             nfilled++;
1259             if ( c1 && c2 ) nempty++;
1260             else c1 = 1;
1261           }
1262           else if ( c1 ) c2 = 1;
1263
1264         }
1265         ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1266       }
1267       ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1268     }
1269   }
1270   return histArray;
1271 }