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