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