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