TPC/AliTPCCalibCE.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 Central Electrode calibration                //
  21 //                                                                                    //
  22 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch  //
  23 //                                                                                    //
  24 ////////////////////////////////////////////////////////////////////////////////////////
  25 //
  26 //
  27 // *************************************************************************************
  28 // *                                Class Description                                  *
  29 // *************************************************************************************
  30 //
  31 /* BEGIN_HTML
  32  <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
  33  using laser runs.</h4>
  34
  35  The information retrieved is
  36  <ul style="list-style-type: square;">
  37    <li>Time arrival from the CE</li>
  38    <li>Signal width</li>
  39    <li>Signal sum</li>
  40  </ul>
  41
  42 <h4>Overview:</h4>
  43  <ol style="list-style-type: upper-roman;">
  44    <li><a href="#working">Working principle</a></li>
  45    <li><a href="#user">User interface for filling data</a></li>
  46    <li><a href="#info">Stored information</a></li>
  47  </ol>
  48
  49  <h3><a name="working">I. Working principle</a></h3>
  50
  51  <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
  52  (see below). These in the end call the Update(...) function.</h4>
  53
  54  <ul style="list-style-type: square;">
  55    <li>the Update(...) function:<br />
  56        In this function the array fPadSignal is filled with the adc signals between the specified range
  57        fFirstTimeBin and fLastTimeBin for the current pad.
  58        before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
  59        stored in fPadSignal.
  60    </li>
  61    <ul style="list-style-type: square;">
  62    <li>the ProcessPad() function:</li>
  63    <ol style="list-style-type: decimal;">
  64      <li>Find Pedestal and Noise information</li>
  65      <ul style="list-style-type: square;">
  66        <li>use database information which has to be set by calling<br />
  67            SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
  68        <li>if no information from the pedestal data base
  69            is available the informaion is calculated on the fly
  70            ( see FindPedestal() function )</li>
  71      </ul>
  72      <li>Find local maxima of the pad signal</li>
  73      <ul style="list-style-type: square;">
  74        <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
  75        have been observed ( see FindLocalMaxima(...) )</li>
  76      </ul>
  77      <li>Find the CE signal information</li>
  78      <ul style="list-style-type: square;">
  79        <li>to find the position of the CE signal the Tmean information from the previos event is used
  80            as the CE signal the local maximum closest to this Tmean is identified</li>
  81        <li>calculate  mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
  82            the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
  83      </ul>
  84      <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
  85      <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
  86      </ol>
  87    </ul>
  88  </ul>
  89
  90  <h4>At the end of each event the EndEvent() function is called</h4>
  91
  92  <ul style="list-style-type: square;">
  93    <li>the EndEvent() function:</li>
  94    <ul style="list-style-type: square;">
  95      <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
  96          This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
  97      <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
  98      <li>calculate Mean Q</li>
  99      <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
 100    </ul>
 101  </ul>
 102
 103  <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
 104   <ul style="list-style-type: square;">
 105   <li>the Analyse() function:</li>
 106     <ul style="list-style-type: square;">
 107       <li>calculate the mean values of T0, RMS, Q for each pad, using
 108           the AliMathBase::GetCOG(...) function</li>
 109       <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
 110          (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
 111     </ul>
 112   </ul>
 113
 114  <h3><a name="user">II. User interface for filling data</a></h3>
 115
 116  <h4>To Fill information one of the following functions can be used:</h4>
 117
 118  <ul style="list-style-type: none;">
 119   <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
 120     <ul style="list-style-type: square;">
 121       <li>process Date event</li>
 122       <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
 123     </ul>
 124     <br />
 125
 126   <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
 127     <ul style="list-style-type: square;">
 128       <li>process AliRawReader event</li>
 129       <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
 130     </ul>
 131     <br />
 132
 133   <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
 134     <ul style="list-style-type: square;">
 135       <li>process event from AliTPCRawStream</li>
 136       <li>call Update function for signal filling</li>
 137     </ul>
 138     <br />
 139
 140   <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 141               iPad,  const Int_t iTimeBin, const Float_t signal);</li>
 142     <ul style="list-style-type: square;">
 143       <li>directly  fill signal information (sector, row, pad, time bin, pad)
 144           to the reference histograms</li>
 145     </ul>
 146  </ul>
 147
 148  <h4>It is also possible to merge two independently taken calibrations using the function</h4>
 149
 150  <ul style="list-style-type: none;">
 151   <li> void Merge(AliTPCCalibSignal *sig)</li>
 152     <ul style="list-style-type: square;">
 153       <li>copy histograms in 'sig' if they do not exist in this instance</li>
 154       <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
 155       <li>After merging call Analyse again!</li>
 156     </ul>
 157  </ul>
 158
 159
 160  <h4>example: filling data using root raw data:</h4>
 161  <pre>
 162  void fillCE(Char_t *filename)
 163  {
 164     rawReader = new AliRawReaderRoot(fileName);
 165     if ( !rawReader ) return;
 166     AliTPCCalibCE *calib = new AliTPCCalibCE;
 167     while (rawReader->NextEvent()){
 168       calib->ProcessEvent(rawReader);
 169     }
 170     calib->Analyse();
 171     calib->DumpToFile("CEData.root");
 172     delete rawReader;
 173     delete calib;
 174  }
 175  </pre>
 176
 177  <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
 178
 179  <h4><a name="info:stored">III.1 Stored information</a></h4>
 180  <ul style="list-style-type: none;">
 181   <li>Histograms:</li>
 182   <ul style="list-style-type: none;">
 183     <li>For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum])
 184         is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
 185         stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
 186   </ul>
 187   <br />
 188
 189  <li>Calibration Data:</li>
 190  <ul style="list-style-type: none;">
 191       <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
 192           the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
 193           fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
 194           is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
 195  </ul>
 196  <br />
 197
 198  <li>For each event the following information is stored:</li>
 199
 200  <ul style="list-style-type: square;">
 201    <li>event time ( TVectorD  fVEventTime )</li>
 202    <li>event id   ( TVectorD  fVEventNumber )</li>
 203    <br />
 204    <li>mean arrival time for each ROC                ( TObjArray fTMeanArrayEvent )</li>
 205    <li>mean Q for each ROC                           ( TObjArray fQMeanArrayEvent )</li>
 206    <li>parameters of a plane fit for each ROC        ( TObjArray fParamArrayEventPol1 )</li>
 207    <li>parameters of a 2D parabola fit for each ROC  ( TObjArray fParamArrayEventPol2 )</li>
 208   </ul>
 209  </ul>
 210
 211  <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
 212  <ul style="list-style-type: none;">
 213   <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
 214     <ul style="list-style-type: square;">
 215       <li>TH2F *GetHistoT0(Int_t sector);</li>
 216       <li>TH2F *GetHistoRMS(Int_t sector);</li>
 217       <li>TH2F *GetHistoQ(Int_t sector);</li>
 218     </ul>
 219     <br />
 220
 221   <li>Accessing the calibration storage objects:</li>
 222     <ul style="list-style-type: square;">
 223       <li>AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values</li>
 224       <li>AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values</li>
 225       <li>AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values</li>
 226     </ul>
 227     <br />
 228
 229   <li>Accessin the event by event information:</li>
 230     <ul style="list-style-type: square;">
 231       <li>The event by event information can be displayed using the</li>
 232       <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
 233       <li>which creates a graph from the specified variables</li>
 234     </ul>
 235   </ul>
 236
 237   <h4>example for visualisation:</h4>
 238   <pre>
 239   //if the file "CEData.root" was created using the above example one could do the following:
 240   TFile fileCE("CEData.root")
 241   AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
 242   ce->GetCalRocT0(0)->Draw("colz");
 243   ce->GetCalRocRMS(0)->Draw("colz");
 244
 245   //or use the AliTPCCalPad functionality:
 246   AliTPCCalPad padT0(ped->GetCalPadT0());
 247   AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
 248   padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
 249   padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
 250
 251   //display event by event information:
 252   //Draw mean arrival time as a function of the event time for oroc sector A00
 253   ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
 254   //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
 255   ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
 256   </pre>
 257 END_HTML */
 258 //////////////////////////////////////////////////////////////////////////////////////
 259
 260
 261 //Root includes
 262 #include <TObjArray.h>
 263 #include <TH1F.h>
 264 #include <TH2S.h>
 265 #include <TString.h>
 266 #include <TVectorF.h>
 267 #include <TVectorD.h>
 268 #include <TMatrixD.h>
 269 #include <TMath.h>
 270 #include <TGraph.h>
 271 #include <TString.h>
 272 #include <TMap.h>
 273 #include <TDirectory.h>
 274 #include <TSystem.h>
 275 #include <TFile.h>
 276
 277 //AliRoot includes
 278 #include "AliLog.h"
 279 #include "AliRawReader.h"
 280 #include "AliRawReaderRoot.h"
 281 #include "AliRawReaderDate.h"
 282 #include "AliRawEventHeaderBase.h"
 283 #include "AliTPCRawStream.h"
 284 #include "AliTPCRawStreamFast.h"
 285 #include "AliTPCcalibDB.h"
 286 #include "AliTPCCalROC.h"
 287 #include "AliTPCCalPad.h"
 288 #include "AliTPCROC.h"
 289 #include "AliTPCParam.h"
 290 #include "AliTPCCalibCE.h"
 291 #include "AliMathBase.h"
 292 #include "TTreeStream.h"
 293
 294 //date
 295 #include "event.h"
 296 ClassImp(AliTPCCalibCE)
 297
 298
 299 AliTPCCalibCE::AliTPCCalibCE() :
 300   AliTPCCalibRawBase(),
 301   fNbinsT0(200),
 302   fXminT0(-5),
 303   fXmaxT0(5),
 304   fNbinsQ(200),
 305   fXminQ(1),
 306   fXmaxQ(40),
 307   fNbinsRMS(100),
 308   fXminRMS(0.1),
 309   fXmaxRMS(5.1),
 310   fPeakDetMinus(2),
 311   fPeakDetPlus(3),
 312   fPeakIntMinus(2),
 313   fPeakIntPlus(2),
 314   fNoiseThresholdMax(5.),
 315   fNoiseThresholdSum(8.),
 316   fIsZeroSuppressed(kFALSE),
 317   fLastSector(-1),
 318   fSecRejectRatio(.4),
 319   fParam(new AliTPCParam),
 320   fPedestalTPC(0x0),
 321   fPadNoiseTPC(0x0),
 322   fPedestalROC(0x0),
 323   fPadNoiseROC(0x0),
 324   fCalRocArrayT0(72),
 325   fCalRocArrayT0Err(72),
 326   fCalRocArrayQ(72),
 327   fCalRocArrayRMS(72),
 328   fCalRocArrayOutliers(72),
 329   fHistoQArray(72),
 330   fHistoT0Array(72),
 331   fHistoRMSArray(72),
 332   fMeanT0rms(0),
 333   fMeanQrms(0),
 334   fMeanRMSrms(0),
 335   fHistoTmean(72),
 336   fParamArrayEventPol1(72),
 337   fParamArrayEventPol2(72),
 338   fTMeanArrayEvent(72),
 339   fQMeanArrayEvent(72),
 340   fVEventTime(1000),
 341   fVEventNumber(1000),
 342   fVTime0SideA(1000),
 343   fVTime0SideC(1000),
 344   fEventId(-1),
 345   fOldRunNumber(0),
 346   fPadTimesArrayEvent(72),
 347   fPadQArrayEvent(72),
 348   fPadRMSArrayEvent(72),
 349   fPadPedestalArrayEvent(72),
 350   fCurrentChannel(-1),
 351   fCurrentSector(-1),
 352   fCurrentRow(-1),
 353   fMaxPadSignal(-1),
 354   fMaxTimeBin(-1),
 355 //   fPadSignal(1024),
 356   fPadPedestal(0),
 357   fPadNoise(0),
 358   fVTime0Offset(72),
 359   fVTime0OffsetCounter(72),
 360   fVMeanQ(72),
 361   fVMeanQCounter(72),
 362   fCurrentCETimeRef(0)
 363 {
 364   //
 365   // AliTPCSignal default constructor
 366   //
 367   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 368   fFirstTimeBin=650;
 369   fLastTimeBin=1000;
 370   fParam->Update();
 371   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 372 }
 373 //_____________________________________________________________________
 374 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 375   AliTPCCalibRawBase(sig),
 376   fNbinsT0(sig.fNbinsT0),
 377   fXminT0(sig.fXminT0),
 378   fXmaxT0(sig.fXmaxT0),
 379   fNbinsQ(sig.fNbinsQ),
 380   fXminQ(sig.fXminQ),
 381   fXmaxQ(sig.fXmaxQ),
 382   fNbinsRMS(sig.fNbinsRMS),
 383   fXminRMS(sig.fXminRMS),
 384   fXmaxRMS(sig.fXmaxRMS),
 385   fPeakDetMinus(sig.fPeakDetMinus),
 386   fPeakDetPlus(sig.fPeakDetPlus),
 387   fPeakIntMinus(sig.fPeakIntMinus),
 388   fPeakIntPlus(sig.fPeakIntPlus),
 389   fNoiseThresholdMax(sig.fNoiseThresholdMax),
 390   fNoiseThresholdSum(sig.fNoiseThresholdSum),
 391   fIsZeroSuppressed(sig.fIsZeroSuppressed),
 392   fLastSector(-1),
 393   fSecRejectRatio(.4),
 394   fParam(new AliTPCParam),
 395   fPedestalTPC(0x0),
 396   fPadNoiseTPC(0x0),
 397   fPedestalROC(0x0),
 398   fPadNoiseROC(0x0),
 399   fCalRocArrayT0(72),
 400   fCalRocArrayT0Err(72),
 401   fCalRocArrayQ(72),
 402   fCalRocArrayRMS(72),
 403   fCalRocArrayOutliers(72),
 404   fHistoQArray(72),
 405   fHistoT0Array(72),
 406   fHistoRMSArray(72),
 407   fMeanT0rms(sig.fMeanT0rms),
 408   fMeanQrms(sig.fMeanQrms),
 409   fMeanRMSrms(sig.fMeanRMSrms),
 410   fHistoTmean(72),
 411   fParamArrayEventPol1(72),
 412   fParamArrayEventPol2(72),
 413   fTMeanArrayEvent(72),
 414   fQMeanArrayEvent(72),
 415   fVEventTime(sig.fVEventTime),
 416   fVEventNumber(sig.fVEventNumber),
 417   fVTime0SideA(sig.fVTime0SideA),
 418   fVTime0SideC(sig.fVTime0SideC),
 419   fEventId(-1),
 420   fOldRunNumber(0),
 421   fPadTimesArrayEvent(72),
 422   fPadQArrayEvent(72),
 423   fPadRMSArrayEvent(72),
 424   fPadPedestalArrayEvent(72),
 425   fCurrentChannel(-1),
 426   fCurrentSector(-1),
 427   fCurrentRow(-1),
 428   fMaxPadSignal(-1),
 429   fMaxTimeBin(-1),
 430 //   fPadSignal(1024),
 431   fPadPedestal(0),
 432   fPadNoise(0),
 433   fVTime0Offset(72),
 434   fVTime0OffsetCounter(72),
 435   fVMeanQ(72),
 436   fVMeanQCounter(72),
 437   fCurrentCETimeRef(0)
 438 {
 439   //
 440   // AliTPCSignal copy constructor
 441   //
 442   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 443
 444   for (Int_t iSec = 0; iSec < 72; ++iSec){
 445     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 446     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 447     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 448     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 449
 450     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 451     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 452     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 453
 454     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 455     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 456     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 457     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 458
 459     if ( hQ != 0x0 ){
 460       TH2S *hNew = new TH2S(*hQ);
 461       hNew->SetDirectory(0);
 462       fHistoQArray.AddAt(hNew,iSec);
 463     }
 464     if ( hT0 != 0x0 ){
 465       TH2S *hNew = new TH2S(*hT0);
 466       hNew->SetDirectory(0);
 467       fHistoT0Array.AddAt(hNew,iSec);
 468     }
 469     if ( hRMS != 0x0 ){
 470       TH2S *hNew = new TH2S(*hRMS);
 471       hNew->SetDirectory(0);
 472       fHistoRMSArray.AddAt(hNew,iSec);
 473     }
 474   }
 475
 476   //copy fit parameters event by event
 477   TObjArray *arr=0x0;
 478   for (Int_t iSec=0; iSec<72; ++iSec){
 479     arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
 480     if ( arr ){
 481       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 482       fParamArrayEventPol1.AddAt(arrEvents, iSec);
 483       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 484         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 485           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 486     }
 487
 488     arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
 489     if ( arr ){
 490       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 491       fParamArrayEventPol2.AddAt(arrEvents, iSec);
 492       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 493         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 494           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 495     }
 496
 497     TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
 498     TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
 499     if ( vMeanTime )
 500       fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
 501     if ( vMeanQ )
 502       fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
 503   }
 504
 505
 506   fVEventTime.ResizeTo(sig.fVEventTime);
 507   fVEventNumber.ResizeTo(sig.fVEventNumber);
 508   fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
 509   fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 510
 511   fParam->Update();
 512 }
 513 //_____________________________________________________________________
 514 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
 515   AliTPCCalibRawBase(),
 516   fNbinsT0(200),
 517   fXminT0(-5),
 518   fXmaxT0(5),
 519   fNbinsQ(200),
 520   fXminQ(1),
 521   fXmaxQ(40),
 522   fNbinsRMS(100),
 523   fXminRMS(0.1),
 524   fXmaxRMS(5.1),
 525   fPeakDetMinus(2),
 526   fPeakDetPlus(3),
 527   fPeakIntMinus(2),
 528   fPeakIntPlus(2),
 529   fNoiseThresholdMax(5.),
 530   fNoiseThresholdSum(8.),
 531   fIsZeroSuppressed(kFALSE),
 532   fLastSector(-1),
 533   fSecRejectRatio(.4),
 534   fParam(new  AliTPCParam),
 535   fPedestalTPC(0x0),
 536   fPadNoiseTPC(0x0),
 537   fPedestalROC(0x0),
 538   fPadNoiseROC(0x0),
 539   fCalRocArrayT0(72),
 540   fCalRocArrayT0Err(72),
 541   fCalRocArrayQ(72),
 542   fCalRocArrayRMS(72),
 543   fCalRocArrayOutliers(72),
 544   fHistoQArray(72),
 545   fHistoT0Array(72),
 546   fHistoRMSArray(72),
 547   fMeanT0rms(0),
 548   fMeanQrms(0),
 549   fMeanRMSrms(0),
 550   fHistoTmean(72),
 551   fParamArrayEventPol1(72),
 552   fParamArrayEventPol2(72),
 553   fTMeanArrayEvent(72),
 554   fQMeanArrayEvent(72),
 555   fVEventTime(1000),
 556   fVEventNumber(1000),
 557   fVTime0SideA(1000),
 558   fVTime0SideC(1000),
 559   fEventId(-1),
 560   fOldRunNumber(0),
 561   fPadTimesArrayEvent(72),
 562   fPadQArrayEvent(72),
 563   fPadRMSArrayEvent(72),
 564   fPadPedestalArrayEvent(72),
 565   fCurrentChannel(-1),
 566   fCurrentSector(-1),
 567   fCurrentRow(-1),
 568   fMaxPadSignal(-1),
 569   fMaxTimeBin(-1),
 570 //   fPadSignal(1024),
 571   fPadPedestal(0),
 572   fPadNoise(0),
 573   fVTime0Offset(72),
 574   fVTime0OffsetCounter(72),
 575   fVMeanQ(72),
 576   fVMeanQCounter(72),
 577   fCurrentCETimeRef(0)
 578 {
 579   //
 580   // constructor which uses a tmap as input to set some specific parameters
 581   //
 582   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 583   fFirstTimeBin=650;
 584   fLastTimeBin=1000;
 585   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 586   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 587   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 588   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 589   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 590   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 591   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 592   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 593   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 594   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 595   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 596   if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
 597   if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
 598   if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
 599   if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
 600   if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
 601   if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
 602   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 603   if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
 604   if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
 605
 606   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 607
 608   fParam->Update();
 609 }
 610
 611 //_____________________________________________________________________
 612 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 613 {
 614   //
 615   // assignment operator
 616   //
 617   if (&source == this) return *this;
 618   new (this) AliTPCCalibCE(source);
 619
 620   return *this;
 621 }
 622 //_____________________________________________________________________
 623 AliTPCCalibCE::~AliTPCCalibCE()
 624 {
 625     //
 626     // destructor
 627     //
 628
 629     fCalRocArrayT0.Delete();
 630     fCalRocArrayT0Err.Delete();
 631     fCalRocArrayQ.Delete();
 632     fCalRocArrayRMS.Delete();
 633     fCalRocArrayOutliers.Delete();
 634
 635     fHistoQArray.Delete();
 636     fHistoT0Array.Delete();
 637     fHistoRMSArray.Delete();
 638
 639     fHistoTmean.Delete();
 640
 641     fParamArrayEventPol1.Delete();
 642     fParamArrayEventPol2.Delete();
 643     fTMeanArrayEvent.Delete();
 644     fQMeanArrayEvent.Delete();
 645
 646     fPadTimesArrayEvent.Delete();
 647     fPadQArrayEvent.Delete();
 648     fPadRMSArrayEvent.Delete();
 649     fPadPedestalArrayEvent.Delete();
 650
 651 //    if ( fHTime0 ) delete fHTime0;
 652     delete fParam;
 653 }
 654 //_____________________________________________________________________
 655 Int_t AliTPCCalibCE::Update(const Int_t icsector,
 656                                 const Int_t icRow,
 657                                 const Int_t icPad,
 658                                 const Int_t icTimeBin,
 659                                 const Float_t csignal)
 660 {
 661   //
 662   // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 663   // no extra analysis necessary. Assumes knowledge of the signal shape!
 664   // assumes that it is looped over consecutive time bins of one pad
 665   //
 666
 667   //temp
 668
 669   if (icRow<0) return 0;
 670   if (icPad<0) return 0;
 671   if (icTimeBin<0) return 0;
 672   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 673
 674   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 675
 676   //init first pad and sector in this event
 677   if ( fCurrentChannel == -1 ) {
 678     fLastSector=-1;
 679     fCurrentChannel = iChannel;
 680     fCurrentSector  = icsector;
 681     fCurrentRow     = icRow;
 682   }
 683
 684   //process last pad if we change to a new one
 685   if ( iChannel != fCurrentChannel ){
 686     ProcessPad();
 687     fLastSector=fCurrentSector;
 688     fCurrentChannel = iChannel;
 689     fCurrentSector  = icsector;
 690     fCurrentRow     = icRow;
 691   }
 692
 693   //fill signals for current pad
 694   fPadSignal[icTimeBin]=csignal;
 695   if ( csignal > fMaxPadSignal ){
 696     fMaxPadSignal = csignal;
 697     fMaxTimeBin   = icTimeBin;
 698   }
 699   return 0;
 700 }
 701 //_____________________________________________________________________
 702 void AliTPCCalibCE::FindPedestal(Float_t part)
 703 {
 704   //
 705     // find pedestal and noise for the current pad. Use either database or
 706     // truncated mean with part*100%
 707   //
 708   Bool_t noPedestal = kTRUE;
 709
 710     //use pedestal database if set
 711   if (fPedestalTPC&&fPadNoiseTPC){
 712         //only load new pedestals if the sector has changed
 713     if ( fCurrentSector!=fLastSector ){
 714       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 715       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 716     }
 717
 718     if ( fPedestalROC&&fPadNoiseROC ){
 719       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
 720       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 721       noPedestal   = kFALSE;
 722     }
 723
 724   }
 725
 726     //if we are not running with pedestal database, or for the current sector there is no information
 727     //available, calculate the pedestal and noise on the fly
 728   if ( noPedestal ) {
 729     fPadPedestal = 0;
 730     fPadNoise    = 0;
 731     if ( fIsZeroSuppressed ) return;
 732     const Int_t kPedMax = 100;  //maximum pedestal value
 733     Float_t  max    =  0;
 734     Float_t  maxPos =  0;
 735     Int_t    median =  -1;
 736     Int_t    count0 =  0;
 737     Int_t    count1 =  0;
 738     //
 739     Float_t padSignal=0;
 740     //
 741     UShort_t histo[kPedMax];
 742     memset(histo,0,kPedMax*sizeof(UShort_t));
 743
 744         //fill pedestal histogram
 745     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 746       padSignal = fPadSignal[i];
 747       if (padSignal<=0) continue;
 748       if (padSignal>max && i>10) {
 749         max = padSignal;
 750         maxPos = i;
 751       }
 752       if (padSignal>kPedMax-1) continue;
 753       histo[int(padSignal+0.5)]++;
 754       count0++;
 755     }
 756         //find median
 757     for (Int_t i=1; i<kPedMax; ++i){
 758       if (count1<count0*0.5) median=i;
 759       count1+=histo[i];
 760     }
 761         // truncated mean
 762     //
 763     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 764     //
 765     for (Int_t idelta=1; idelta<10; ++idelta){
 766       if (median-idelta<=0) continue;
 767       if (median+idelta>kPedMax) continue;
 768       if (count<part*count1){
 769         count+=histo[median-idelta];
 770         mean +=histo[median-idelta]*(median-idelta);
 771         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 772         count+=histo[median+idelta];
 773         mean +=histo[median+idelta]*(median+idelta);
 774         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 775       }
 776     }
 777     if ( count > 0 ) {
 778       mean/=count;
 779       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 780       fPadPedestal = mean;
 781       fPadNoise    = rms;
 782     }
 783   }
 784 }
 785 //_____________________________________________________________________
 786 void AliTPCCalibCE::UpdateCETimeRef()
 787 {
 788   // Find the time reference of the last valid CE signal in sector
 789   // for irocs of the A-Side the reference of the corresponging OROC is returned
 790   // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
 791   if ( fLastSector == fCurrentSector ) return;
 792   Int_t sector=fCurrentSector;
 793   if ( sector < 18 ) sector+=36;
 794   fCurrentCETimeRef=0;
 795   TVectorF *vtRef = GetTMeanEvents(sector);
 796   if ( !vtRef ) return;
 797   Int_t vtRefSize= vtRef->GetNrows();
 798   if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
 799   else vtRefSize=fNevents;
 800   while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
 801   fCurrentCETimeRef=(*vtRef)[vtRefSize];
 802   AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
 803 }
 804 //_____________________________________________________________________
 805 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
 806 {
 807   //
 808     //  Find position, signal width and height of the CE signal (last signal)
 809     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 810     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 811   //
 812
 813   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 814   Int_t   cemaxpos       = 0;
 815   Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
 816   const Int_t    kCemin  = fPeakIntMinus;             // range for the analysis of the ce signal +- channels from the peak
 817   const Int_t    kCemax  = fPeakIntPlus;
 818
 819   Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
 820
 821     // find maximum closest to the sector mean from the last event
 822   for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
 823         // get sector mean of last event
 824     Float_t tmean = fCurrentCETimeRef;
 825     if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
 826       minDist  = tmean-maxima[imax];
 827       cemaxpos = (Int_t)maxima[imax];
 828     }
 829   }
 830 //   printf("L1 phase TB: %f\n",GetL1PhaseTB());
 831   if (cemaxpos!=0){
 832     ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
 833     for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
 834       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
 835         Float_t signal = fPadSignal[i]-fPadPedestal;
 836         if (signal>0) {
 837           ceTime+=signal*(i+0.5);
 838           ceRMS +=signal*(i+0.5)*(i+0.5);
 839           ceQsum+=signal;
 840         }
 841       }
 842     }
 843   }
 844   if (ceQmax&&ceQsum>ceSumThreshold) {
 845     ceTime/=ceQsum;
 846     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 847     ceTime-=GetL1PhaseTB();
 848     fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 849     fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 850
 851   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
 852   //                                the pick-up signal should scale with the pad area. In addition
 853   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
 854   //                                ratio 2/3. The pad area we express in cm2. We normalise the signal
 855   //                                to the OROC signal (factor 2/3 for the IROCs).
 856     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 857     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
 858
 859     ceQsum/=norm;
 860     fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
 861     fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
 862   } else {
 863     ceQmax=0;
 864     ceTime=0;
 865     ceRMS =0;
 866     ceQsum=0;
 867   }
 868   param[0] = ceQmax;
 869   param[1] = ceTime;
 870   param[2] = ceRMS;
 871   qSum     = ceQsum;
 872 }
 873 //_____________________________________________________________________
 874 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
 875 {
 876   //
 877   // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
 878   // and 'tplus' timebins after 'pos'
 879   //
 880   if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
 881   for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
 882     if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
 883   for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
 884     if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
 885       ++iTime2;
 886     if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
 887   }
 888   return kTRUE;
 889 }
 890 //_____________________________________________________________________
 891 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
 892 {
 893   //
 894     // Find local maxima on the pad signal and Histogram them
 895   //
 896   Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
 897   Int_t   count       = 0;
 898
 899   for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
 900     if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
 901     if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
 902       if (count<maxima.GetNrows()){
 903         maxima.GetMatrixArray()[count++]=i;
 904         GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
 905         i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin  fPeakDetMinus+fPeakDetPlus-1
 906       }
 907     }
 908   }
 909 }
 910 //_____________________________________________________________________
 911 void AliTPCCalibCE::ProcessPad()
 912 {
 913   //
 914   //  Process data of current pad
 915   //
 916   FindPedestal();
 917
 918   TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
 919                              // + central electrode and possibly post peaks from the CE signal
 920                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
 921   FindLocalMaxima(maxima);
 922   if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
 923
 924   UpdateCETimeRef();                       // update the time refenrence for the current sector
 925   if ( fCurrentCETimeRef==0 ) return;      //return if we don't have time 0 info, eg if only one side has laser
 926   TVectorD param(3);
 927   Float_t  qSum;
 928   FindCESignal(param, qSum, maxima);
 929
 930   Double_t meanT  = param[1];
 931   Double_t sigmaT = param[2];
 932
 933     //Fill Event T0 counter
 934   (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 935
 936     //Fill Q histogram
 937   GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
 938
 939     //Fill RMS histogram
 940   GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 941
 942
 943     //Fill debugging info
 944   if ( GetStreamLevel()>0 ){
 945     (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
 946     (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
 947     (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
 948   }
 949
 950   ResetPad();
 951 }
 952 //_____________________________________________________________________
 953 void AliTPCCalibCE::EndEvent()
 954 {
 955   //  Process data of current pad
 956   //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
 957   //  before the EndEvent function to set the event timestamp and number!!!
 958   //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
 959   //  function was called
 960
 961   //check if last pad has allready been processed, if not do so
 962   if ( fMaxTimeBin>-1 ) ProcessPad();
 963
 964   AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
 965
 966   TVectorD param(3);
 967   TMatrixD dummy(3,3);
 968 //    TVectorF vMeanTime(72);
 969 //    TVectorF vMeanQ(72);
 970   AliTPCCalROC *calIroc=new AliTPCCalROC(0);
 971   AliTPCCalROC *calOroc=new AliTPCCalROC(36);
 972
 973   //find mean time0 offset for side A and C
 974   //use only orocs due to the better statistics
 975   Double_t time0Side[2];       //time0 for side A:0 and C:1
 976   Double_t time0SideCount[2];  //time0 counter for side A:0 and C:1
 977   time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
 978   for ( Int_t iSec = 36; iSec<72; ++iSec ){
 979     time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
 980     time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
 981   }
 982   if ( time0SideCount[0] >0  )
 983     time0Side[0]/=time0SideCount[0];
 984   if ( time0SideCount[1] >0 )
 985     time0Side[1]/=time0SideCount[1];
 986     // end find time0 offset
 987   AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
 988   Int_t nSecMeanT=0;
 989   //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
 990   for ( Int_t iSec = 0; iSec<72; ++iSec ){
 991     AliDebug(4,Form("Processing sector '%02d'\n",iSec));
 992     //find median and then calculate the mean around it
 993     TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
 994     if ( !hMeanT ) continue;
 995     //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
 996     if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
 997       hMeanT->Reset();
 998       AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
 999       continue;
1000     }
1001
1002     Double_t entries = hMeanT->GetEffectiveEntries();
1003     Double_t sum     = 0;
1004     Short_t *arr     = hMeanT->GetArray()+1;
1005     Int_t ibin=0;
1006     for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1007       sum+=arr[ibin];
1008       if ( sum>=(entries/2.) ) break;
1009     }
1010     Int_t delta = 4;
1011     Int_t firstBin = fFirstTimeBin+ibin-delta;
1012     Int_t lastBin  = fFirstTimeBin+ibin+delta;
1013     if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1014     if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
1015     Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1016
1017         // check boundaries for ebye info of mean time
1018     TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1019     Int_t vSize=vMeanTime->GetNrows();
1020     if ( vSize < fNevents+1 ){
1021       vMeanTime->ResizeTo(vSize+100);
1022     }
1023
1024     // store mean time for the readout sides
1025     vSize=fVTime0SideA.GetNrows();
1026     if ( vSize < fNevents+1 ){
1027       fVTime0SideA.ResizeTo(vSize+100);
1028       fVTime0SideC.ResizeTo(vSize+100);
1029     }
1030     fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1031     fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1032
1033     vMeanTime->GetMatrixArray()[fNevents]=median;
1034     nSecMeanT++;
1035     // end find median
1036
1037     TVectorF *vTimes = GetPadTimesEvent(iSec);
1038     if ( !vTimes ) continue;                     //continue if no time information for this sector is available
1039
1040     AliTPCCalROC calIrocOutliers(0);
1041     AliTPCCalROC calOrocOutliers(36);
1042
1043     // calculate mean Q of the sector
1044     TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1045     vSize=vMeanQ->GetNrows();
1046     if ( vSize < fNevents+1 ){
1047       vMeanQ->ResizeTo(vSize+100);
1048     }
1049     Float_t meanQ = 0;
1050     if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1051     vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1052
1053     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1054       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
1055
1056             //set values for temporary roc calibration class
1057       if ( iSec < 36 ) {
1058         calIroc->SetValue(iChannel, time);
1059         if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
1060
1061       } else {
1062         calOroc->SetValue(iChannel, time);
1063         if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
1064       }
1065
1066       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1067         // test that we really found the CE signal reliably
1068         if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1069           GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1070
1071
1072
1073             //-------------------------------  Debug start  ------------------------------
1074       if ( GetStreamLevel()>0 ){
1075         TTreeSRedirector *streamer=GetDebugStreamer();
1076         if (streamer){
1077           Int_t row=0;
1078           Int_t pad=0;
1079           Int_t padc=0;
1080
1081           Float_t q   = (*GetPadQEvent(iSec))[iChannel];
1082           Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1083
1084           UInt_t channel=iChannel;
1085           Int_t sector=iSec;
1086
1087           while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1088           pad = channel-fROC->GetRowIndexes(sector)[row];
1089           padc = pad-(fROC->GetNPads(sector,row)/2);
1090
1091 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1092 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
1093 //                                  fLastTimeBin-fFirstTimeBin,
1094 //                                  fFirstTimeBin,fLastTimeBin);
1095 //              h1->SetDirectory(0);
1096         //
1097 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1098 //                  h1->Fill(i,fPadSignal(i));
1099
1100           Double_t t0Sec = 0;
1101           if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1102             t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1103           Double_t t0Side = time0Side[(iSec/18)%2];
1104           (*streamer) << "DataPad" <<
1105             "Event=" << fNevents <<
1106             "RunNumber=" << fRunNumber <<
1107             "TimeStamp="   << fTimeStamp <<
1108             "Sector="<< sector <<
1109             "Row="   << row<<
1110             "Pad="   << pad <<
1111             "PadC="  << padc <<
1112             "PadSec="<< channel <<
1113             "Time0Sec="  << t0Sec <<
1114             "Time0Side=" << t0Side <<
1115             "Time="  << time <<
1116             "RMS="   << rms <<
1117             "Sum="   << q <<
1118             "MeanQ=" << meanQ <<
1119         //                  "hist.=" << h1 <<
1120             "\n";
1121
1122     //          delete h1;
1123         }
1124       }
1125       //-----------------------------  Debug end  ------------------------------
1126     }// end channel loop
1127
1128
1129     //do fitting now only in debug mode
1130     if (GetDebugLevel()>0){
1131       TVectorD paramPol1(3);
1132       TVectorD paramPol2(6);
1133       TMatrixD matPol1(3,3);
1134       TMatrixD matPol2(6,6);
1135       Float_t  chi2Pol1=0;
1136       Float_t  chi2Pol2=0;
1137
1138       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1139         if ( iSec < 36 ){
1140           calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1141           calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1142         } else {
1143           calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1144           calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1145         }
1146
1147         GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1148         GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1149       }
1150
1151   //-------------------------------  Debug start  ------------------------------
1152       if ( GetStreamLevel()>0 ){
1153         TTreeSRedirector *streamer=GetDebugStreamer();
1154         if ( streamer ) {
1155           (*streamer) << "DataRoc" <<
1156 //    "Event=" << fEvent <<
1157             "RunNumber=" << fRunNumber <<
1158             "TimeStamp="   << fTimeStamp <<
1159             "Sector="<< iSec <<
1160             "hMeanT.=" << hMeanT <<
1161             "median=" << median <<
1162             "paramPol1.=" << &paramPol1 <<
1163             "paramPol2.=" << &paramPol2 <<
1164             "matPol1.="   << &matPol1 <<
1165             "matPol2.="   << &matPol2 <<
1166             "chi2Pol1="   << chi2Pol1 <<
1167             "chi2Pol2="   << chi2Pol2 <<
1168             "\n";
1169         }
1170       }
1171     }
1172         //-------------------------------  Debug end  ------------------------------
1173     hMeanT->Reset();
1174   }// end sector loop
1175     //return if no sector has a valid mean time
1176   if ( nSecMeanT == 0 ) return;
1177
1178
1179 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1180 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1181   if ( fVEventTime.GetNrows() < fNevents+1 ) {
1182     fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1183     fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1184   }
1185   fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1186   fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1187
1188   fNevents++;
1189   fOldRunNumber = fRunNumber;
1190
1191   delete calIroc;
1192   delete calOroc;
1193   AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1194 }
1195 //_____________________________________________________________________
1196 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1197                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1198                                   const Char_t *type, Bool_t force)
1199 {
1200     //
1201     // return pointer to TH2S histogram of 'type'
1202     // if force is true create a new histogram if it doesn't exist allready
1203     //
1204     if ( !force || arr->UncheckedAt(sector) )
1205         return (TH2S*)arr->UncheckedAt(sector);
1206
1207     // if we are forced and histogram doesn't exist yet create it
1208     Char_t name[255], title[255];
1209
1210     sprintf(name,"hCalib%s%.2d",type,sector);
1211     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1212
1213     // new histogram with Q calib information. One value for each pad!
1214     TH2S* hist = new TH2S(name,title,
1215                           nbinsY, ymin, ymax,
1216                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1217     hist->SetDirectory(0);
1218     arr->AddAt(hist,sector);
1219     return hist;
1220 }
1221 //_____________________________________________________________________
1222 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1223 {
1224     //
1225     // return pointer to T0 histogram
1226     // if force is true create a new histogram if it doesn't exist allready
1227     //
1228     TObjArray *arr = &fHistoT0Array;
1229     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1230 }
1231 //_____________________________________________________________________
1232 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1233 {
1234     //
1235     // return pointer to Q histogram
1236     // if force is true create a new histogram if it doesn't exist allready
1237     //
1238     TObjArray *arr = &fHistoQArray;
1239     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1240 }
1241 //_____________________________________________________________________
1242 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1243 {
1244     //
1245     // return pointer to Q histogram
1246     // if force is true create a new histogram if it doesn't exist allready
1247     //
1248     TObjArray *arr = &fHistoRMSArray;
1249     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1250 }
1251 //_____________________________________________________________________
1252 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1253                               const Char_t *type, Bool_t force)
1254 {
1255     //
1256     // return pointer to TH1S histogram
1257     // if force is true create a new histogram if it doesn't exist allready
1258     //
1259     if ( !force || arr->UncheckedAt(sector) )
1260         return (TH1S*)arr->UncheckedAt(sector);
1261
1262     // if we are forced and histogram doesn't yes exist create it
1263     Char_t name[255], title[255];
1264
1265     sprintf(name,"hCalib%s%.2d",type,sector);
1266     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1267
1268     // new histogram with calib information. One value for each pad!
1269     TH1S* hist = new TH1S(name,title,
1270                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1271     hist->SetDirectory(0);
1272     arr->AddAt(hist,sector);
1273     return hist;
1274 }
1275 //_____________________________________________________________________
1276 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1277 {
1278     //
1279     // return pointer to Q histogram
1280     // if force is true create a new histogram if it doesn't exist allready
1281     //
1282     TObjArray *arr = &fHistoTmean;
1283     return GetHisto(sector, arr, "LastTmean", force);
1284 }
1285 //_____________________________________________________________________
1286 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1287 {
1288     //
1289     // return pointer to Pad Info from 'arr' for the current event and sector
1290     // if force is true create it if it doesn't exist allready
1291     //
1292     if ( !force || arr->UncheckedAt(sector) )
1293         return (TVectorF*)arr->UncheckedAt(sector);
1294
1295     TVectorF *vect = new TVectorF(size);
1296     arr->AddAt(vect,sector);
1297     return vect;
1298 }
1299 //_____________________________________________________________________
1300 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1301 {
1302     //
1303     // return pointer to Pad Times Array for the current event and sector
1304     // if force is true create it if it doesn't exist allready
1305     //
1306     TObjArray *arr = &fPadTimesArrayEvent;
1307     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1308 }
1309 //_____________________________________________________________________
1310 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1311 {
1312     //
1313     // return pointer to Pad Q Array for the current event and sector
1314     // if force is true create it if it doesn't exist allready
1315     // for debugging purposes only
1316     //
1317
1318     TObjArray *arr = &fPadQArrayEvent;
1319     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1320 }
1321 //_____________________________________________________________________
1322 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1323 {
1324     //
1325     // return pointer to Pad RMS Array for the current event and sector
1326     // if force is true create it if it doesn't exist allready
1327     // for debugging purposes only
1328     //
1329     TObjArray *arr = &fPadRMSArrayEvent;
1330     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1331 }
1332 //_____________________________________________________________________
1333 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1334 {
1335     //
1336     // return pointer to Pad RMS Array for the current event and sector
1337     // if force is true create it if it doesn't exist allready
1338     // for debugging purposes only
1339     //
1340     TObjArray *arr = &fPadPedestalArrayEvent;
1341     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1342 }
1343 //_____________________________________________________________________
1344 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1345 {
1346     //
1347     // return pointer to the EbyE info of the mean arrival time for 'sector'
1348     // if force is true create it if it doesn't exist allready
1349     //
1350     TObjArray *arr = &fTMeanArrayEvent;
1351     return GetVectSector(sector,arr,100,force);
1352 }
1353 //_____________________________________________________________________
1354 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1355 {
1356     //
1357     // return pointer to the EbyE info of the mean arrival time for 'sector'
1358     // if force is true create it if it doesn't exist allready
1359     //
1360     TObjArray *arr = &fQMeanArrayEvent;
1361     return GetVectSector(sector,arr,100,force);
1362 }
1363 //_____________________________________________________________________
1364 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1365 {
1366     //
1367     // return pointer to ROC Calibration
1368     // if force is true create a new histogram if it doesn't exist allready
1369     //
1370     if ( !force || arr->UncheckedAt(sector) )
1371         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1372
1373     // if we are forced and histogram doesn't yes exist create it
1374
1375     // new AliTPCCalROC for T0 information. One value for each pad!
1376     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1377     arr->AddAt(croc,sector);
1378     return croc;
1379 }
1380 //_____________________________________________________________________
1381 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1382 {
1383     //
1384     // return pointer to Time 0 ROC Calibration
1385     // if force is true create a new histogram if it doesn't exist allready
1386     //
1387     TObjArray *arr = &fCalRocArrayT0;
1388     return GetCalRoc(sector, arr, force);
1389 }
1390 //_____________________________________________________________________
1391 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1392 {
1393     //
1394     // return pointer to the error of Time 0 ROC Calibration
1395     // if force is true create a new histogram if it doesn't exist allready
1396     //
1397     TObjArray *arr = &fCalRocArrayT0Err;
1398     return GetCalRoc(sector, arr, force);
1399 }
1400 //_____________________________________________________________________
1401 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1402 {
1403     //
1404     // return pointer to T0 ROC Calibration
1405     // if force is true create a new histogram if it doesn't exist allready
1406     //
1407     TObjArray *arr = &fCalRocArrayQ;
1408     return GetCalRoc(sector, arr, force);
1409 }
1410 //_____________________________________________________________________
1411 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1412 {
1413     //
1414     // return pointer to signal width ROC Calibration
1415     // if force is true create a new histogram if it doesn't exist allready
1416     //
1417     TObjArray *arr = &fCalRocArrayRMS;
1418     return GetCalRoc(sector, arr, force);
1419 }
1420 //_____________________________________________________________________
1421 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1422 {
1423     //
1424     // return pointer to Outliers
1425     // if force is true create a new histogram if it doesn't exist allready
1426     //
1427     TObjArray *arr = &fCalRocArrayOutliers;
1428     return GetCalRoc(sector, arr, force);
1429 }
1430 //_____________________________________________________________________
1431 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1432 {
1433     //
1434     // return pointer to TObjArray of fit parameters
1435     // if force is true create a new histogram if it doesn't exist allready
1436     //
1437     if ( !force || arr->UncheckedAt(sector) )
1438         return (TObjArray*)arr->UncheckedAt(sector);
1439
1440     // if we are forced and array doesn't yes exist create it
1441
1442     // new TObjArray for parameters
1443     TObjArray *newArr = new TObjArray;
1444     arr->AddAt(newArr,sector);
1445     return newArr;
1446 }
1447 //_____________________________________________________________________
1448 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1449 {
1450     //
1451     // return pointer to TObjArray of fit parameters from plane fit
1452     // if force is true create a new histogram if it doesn't exist allready
1453     //
1454     TObjArray *arr = &fParamArrayEventPol1;
1455     return GetParamArray(sector, arr, force);
1456 }
1457 //_____________________________________________________________________
1458 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1459 {
1460     //
1461     // return pointer to TObjArray of fit parameters from parabola fit
1462     // if force is true create a new histogram if it doesn't exist allready
1463     //
1464     TObjArray *arr = &fParamArrayEventPol2;
1465     return GetParamArray(sector, arr, force);
1466 }
1467 //_____________________________________________________________________
1468 void AliTPCCalibCE::ResetEvent()
1469 {
1470     //
1471     //  Reset global counters  -- Should be called before each event is processed
1472     //
1473     fLastSector=-1;
1474     fCurrentSector=-1;
1475     fCurrentRow=-1;
1476     fCurrentChannel=-1;
1477
1478     ResetPad();
1479
1480     fPadTimesArrayEvent.Delete();
1481     fPadQArrayEvent.Delete();
1482     fPadRMSArrayEvent.Delete();
1483     fPadPedestalArrayEvent.Delete();
1484
1485     for ( Int_t i=0; i<72; ++i ){
1486         fVTime0Offset.GetMatrixArray()[i]=0;
1487         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1488         fVMeanQ.GetMatrixArray()[i]=0;
1489         fVMeanQCounter.GetMatrixArray()[i]=0;
1490     }
1491 }
1492 //_____________________________________________________________________
1493 void AliTPCCalibCE::ResetPad()
1494 {
1495     //
1496     //  Reset pad infos -- Should be called after a pad has been processed
1497     //
1498     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1499       fPadSignal[i] = 0;
1500     fMaxTimeBin   = -1;
1501     fMaxPadSignal = -1;
1502     fPadPedestal  = -1;
1503     fPadNoise     = -1;
1504 }
1505 //_____________________________________________________________________
1506 void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1507 {
1508     //
1509     //  Merge ce to the current AliTPCCalibCE
1510     //
1511
1512     //merge histograms
1513     for (Int_t iSec=0; iSec<72; ++iSec){
1514         TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1515         TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1516         TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1517
1518
1519         if ( hRefQmerge ){
1520             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1521             TH2S *hRefQ   = GetHistoQ(iSec);
1522             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1523             else {
1524                 TH2S *hist = new TH2S(*hRefQmerge);
1525                 hist->SetDirectory(0);
1526                 fHistoQArray.AddAt(hist, iSec);
1527             }
1528             hRefQmerge->SetDirectory(dir);
1529         }
1530         if ( hRefT0merge ){
1531             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1532             TH2S *hRefT0  = GetHistoT0(iSec);
1533             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1534             else {
1535                 TH2S *hist = new TH2S(*hRefT0merge);
1536                 hist->SetDirectory(0);
1537                 fHistoT0Array.AddAt(hist, iSec);
1538             }
1539             hRefT0merge->SetDirectory(dir);
1540         }
1541         if ( hRefRMSmerge ){
1542             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1543             TH2S *hRefRMS = GetHistoRMS(iSec);
1544             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1545             else {
1546                 TH2S *hist = new TH2S(*hRefRMSmerge);
1547                 hist->SetDirectory(0);
1548                 fHistoRMSArray.AddAt(hist, iSec);
1549             }
1550             hRefRMSmerge->SetDirectory(dir);
1551         }
1552
1553     }
1554
1555     // merge time information
1556
1557
1558     Int_t nCEevents = ce->GetNeventsProcessed();
1559     for (Int_t iSec=0; iSec<72; ++iSec){
1560         TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1561         TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1562         TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1563         TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1564
1565         TObjArray *arrPol1  = 0x0;
1566         TObjArray *arrPol2  = 0x0;
1567         TVectorF *vMeanTime = 0x0;
1568         TVectorF *vMeanQ    = 0x0;
1569
1570         //resize arrays
1571         if ( arrPol1CE && arrPol2CE ){
1572             arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1573             arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1574             arrPol1->Expand(fNevents+nCEevents);
1575             arrPol2->Expand(fNevents+nCEevents);
1576         }
1577         if ( vMeanTimeCE && vMeanQCE ){
1578             vMeanTime = GetTMeanEvents(iSec,kTRUE);
1579             vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1580             vMeanTime->ResizeTo(fNevents+nCEevents);
1581             vMeanQ->ResizeTo(fNevents+nCEevents);
1582         }
1583
1584
1585         for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1586             if ( arrPol1CE && arrPol2CE ){
1587                 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1588                 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1589                 if ( paramPol1 && paramPol2 ){
1590                     GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1591                     GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1592                 }
1593             }
1594             if ( vMeanTimeCE && vMeanQCE ){
1595                 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1596                 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1597             }
1598         }
1599     }
1600
1601
1602
1603     TVectorD*  eventTimes  = ce->GetEventTimes();
1604     TVectorD*  eventIds  = ce->GetEventIds();
1605     fVEventTime.ResizeTo(fNevents+nCEevents);
1606     fVEventNumber.ResizeTo(fNevents+nCEevents);
1607
1608     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1609         Double_t evTime     = eventTimes->GetMatrixArray()[iEvent];
1610         Double_t evId       = eventIds->GetMatrixArray()[iEvent];
1611         fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1612         fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1613     }
1614     fNevents+=nCEevents; //increase event counter
1615
1616 }
1617 //_____________________________________________________________________
1618 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1619 {
1620   //
1621   // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1622   // or side (-1: A-Side, -2: C-Side)
1623   // xVariable:    0-event time, 1-event id, 2-internal event counter
1624   // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1625   // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1626   //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1627   //                               not used for mean time and mean Q )
1628   // for an example see class description at the beginning
1629   //
1630
1631   Double_t *x = new Double_t[fNevents];
1632   Double_t *y = new Double_t[fNevents];
1633
1634   TVectorD *xVar = 0x0;
1635   TObjArray *aType = 0x0;
1636   Int_t npoints=0;
1637
1638     // sanity checks
1639   if ( (sector<-2) || (sector>71) )      return 0x0;
1640   if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1641   if ( (fitType<0) || (fitType>3) )     return 0x0;
1642   if ( sector>=0&&!GetTMeanEvents(sector) )        return 0x0; //no mean time information available
1643   if ( sector<0 && fitType!=2) return 0x0;
1644
1645   if (sector>=0){
1646     if ( fitType==0 ){
1647       if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1648       aType = &fParamArrayEventPol1;
1649       if ( aType->At(sector)==0x0 ) return 0x0;
1650     }
1651     else if ( fitType==1 ){
1652       if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1653       aType = &fParamArrayEventPol2;
1654       if ( aType->At(sector)==0x0 ) return 0x0;
1655     }
1656
1657   }
1658   if ( xVariable == 0 ) xVar = &fVEventTime;
1659   if ( xVariable == 1 ) xVar = &fVEventNumber;
1660   if ( xVariable == 2 ) {
1661     xVar = new TVectorD(fNevents);
1662     for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1663   }
1664
1665   for (Int_t ievent =0; ievent<fNevents; ++ievent){
1666     if ( fitType<2 ){
1667       TObjArray *events = (TObjArray*)(aType->At(sector));
1668       if ( events->GetSize()<=ievent ) break;
1669       TVectorD *v = (TVectorD*)(events->At(ievent));
1670       if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1671     } else if (fitType == 2) {
1672       Double_t xValue=(*xVar)[ievent];
1673       Double_t yValue=0;
1674       if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1675       else if (sector==-1) yValue=fVTime0SideA(ievent);
1676       else if (sector==-2) yValue=fVTime0SideC(ievent);
1677       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1678     }else if (fitType == 3) {
1679       Double_t xValue=(*xVar)[ievent];
1680       Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1681       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1682     }
1683   }
1684
1685   TGraph *gr = new TGraph(npoints);
1686     //sort xVariable increasing
1687   Int_t    *sortIndex = new Int_t[npoints];
1688   TMath::Sort(npoints,x,sortIndex);
1689   for (Int_t i=0;i<npoints;++i){
1690     gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1691   }
1692
1693
1694   if ( xVariable == 2 ) delete xVar;
1695   delete x;
1696   delete y;
1697   delete sortIndex;
1698   return gr;
1699 }
1700 //_____________________________________________________________________
1701 void AliTPCCalibCE::Analyse()
1702 {
1703   //
1704   //  Calculate calibration constants
1705   //
1706
1707   TVectorD paramQ(3);
1708   TVectorD paramT0(3);
1709   TVectorD paramRMS(3);
1710   TMatrixD dummy(3,3);
1711
1712   Float_t channelCounter=0;
1713   fMeanT0rms=0;
1714   fMeanQrms=0;
1715   fMeanRMSrms=0;
1716
1717   for (Int_t iSec=0; iSec<72; ++iSec){
1718     TH2S *hT0 = GetHistoT0(iSec);
1719     if (!hT0 ) continue;
1720
1721     AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1722     AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1723     AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1724     AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1725     AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
1726
1727     TH2S *hQ   = GetHistoQ(iSec);
1728     TH2S *hRMS = GetHistoRMS(iSec);
1729
1730     Short_t *arrayhQ   = hQ->GetArray();
1731     Short_t *arrayhT0  = hT0->GetArray();
1732     Short_t *arrayhRMS = hRMS->GetArray();
1733
1734     UInt_t nChannels = fROC->GetNChannels(iSec);
1735
1736   //debug
1737     Int_t row=0;
1738     Int_t pad=0;
1739     Int_t padc=0;
1740   //! debug
1741
1742     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1743
1744
1745       Float_t cogTime0 = -1000;
1746       Float_t cogQ     = -1000;
1747       Float_t cogRMS   = -1000;
1748       Float_t cogOut   = 0;
1749       Float_t rms      = 0;
1750       Float_t rmsT0    = 0;
1751
1752
1753       Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1754       Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1755       Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1756
1757       cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
1758       fMeanQrms+=rms;
1759       cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
1760       fMeanT0rms+=rmsT0;
1761       cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
1762       fMeanRMSrms+=rms;
1763       channelCounter++;
1764
1765       /*
1766              //outlier specifications
1767       if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1768     cogOut = 1;
1769     cogTime0 = 0;
1770     cogQ     = 0;
1771     cogRMS   = 0;
1772       }
1773 */
1774       rocQ->SetValue(iChannel, cogQ*cogQ);
1775       rocT0->SetValue(iChannel, cogTime0);
1776       rocT0Err->SetValue(iChannel, rmsT0);
1777       rocRMS->SetValue(iChannel, cogRMS);
1778       rocOut->SetValue(iChannel, cogOut);
1779
1780
1781       //debug
1782       if ( GetStreamLevel() > 0 ){
1783         TTreeSRedirector *streamer=GetDebugStreamer();
1784         if ( streamer ) {
1785
1786           while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1787           pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1788           padc = pad-(fROC->GetNPads(iSec,row)/2);
1789
1790           (*streamer) << "DataEnd" <<
1791             "Sector="  << iSec      <<
1792             "Pad="     << pad       <<
1793             "PadC="    << padc      <<
1794             "Row="     << row       <<
1795             "PadSec="  << iChannel   <<
1796             "Q="       << cogQ      <<
1797             "T0="      << cogTime0  <<
1798             "RMS="     << cogRMS    <<
1799             "\n";
1800         }
1801       }
1802       //! debug
1803
1804     }
1805
1806   }
1807   if ( channelCounter>0 ){
1808     fMeanT0rms/=channelCounter;
1809     fMeanQrms/=channelCounter;
1810     fMeanRMSrms/=channelCounter;
1811   }
1812 //   if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1813 //    delete fDebugStreamer;
1814 //    fDebugStreamer = 0x0;
1815   fVEventTime.ResizeTo(fNevents);
1816   fVEventNumber.ResizeTo(fNevents);
1817   fVTime0SideA.ResizeTo(fNevents);
1818   fVTime0SideC.ResizeTo(fNevents);
1819 }