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