TPC/AliTPCCalibPulser.cxx

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