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