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