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