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 <TH1.h>
 264 #include <TH1F.h>
 265 #include <TH2S.h>
 266 #include <TF1.h>
 267 #include <TString.h>
 268 #include <TVectorF.h>
 269 #include <TVectorD.h>
 270 #include <TVector3.h>
 271 #include <TMatrixD.h>
 272 #include <TMath.h>
 273 #include <TGraph.h>
 274 #include <TGraphErrors.h>
 275 #include <TString.h>
 276 #include <TMap.h>
 277 #include <TDirectory.h>
 278 #include <TSystem.h>
 279 #include <TFile.h>
 280 #include <TCollection.h>
 281 #include <TTimeStamp.h>
 282 #include <TList.h>
 283 #include <TKey.h>
 284 #include <TSpectrum.h>
 285
 286 //AliRoot includes
 287 #include "AliLog.h"
 288 #include "AliRawReader.h"
 289 #include "AliRawReaderRoot.h"
 290 #include "AliRawReaderDate.h"
 291 #include "AliRawEventHeaderBase.h"
 292 #include "AliTPCRawStream.h"
 293 #include "AliTPCCalROC.h"
 294 #include "AliTPCCalPad.h"
 295 #include "AliTPCROC.h"
 296 #include "AliTPCParam.h"
 297 #include "AliTPCCalibCE.h"
 298 #include "AliMathBase.h"
 299 #include "AliTPCTransform.h"
 300 #include "AliTPCLaserTrack.h"
 301 #include "TTreeStream.h"
 302
 303 #include "AliCDBManager.h"
 304 #include "AliCDBEntry.h"
 305 //date
 306 #include "event.h"
 307 ClassImp(AliTPCCalibCE)
 308
 309
 310 AliTPCCalibCE::AliTPCCalibCE() :
 311   AliTPCCalibRawBase(),
 312   fNbinsT0(200),
 313   fXminT0(-5),
 314   fXmaxT0(5),
 315   fNbinsQ(200),
 316   fXminQ(1),
 317   fXmaxQ(40),
 318   fNbinsRMS(100),
 319   fXminRMS(0.1),
 320   fXmaxRMS(5.1),
 321   fPeakDetMinus(2),
 322   fPeakDetPlus(3),
 323   fPeakIntMinus(2),
 324   fPeakIntPlus(2),
 325   fNoiseThresholdMax(5.),
 326   fNoiseThresholdSum(8.),
 327   fIsZeroSuppressed(kFALSE),
 328   fLastSector(-1),
 329   fSecRejectRatio(.4),
 330   fParam(new AliTPCParam),
 331   fPedestalTPC(0x0),
 332   fPadNoiseTPC(0x0),
 333   fPedestalROC(0x0),
 334   fPadNoiseROC(0x0),
 335   fCalRocArrayT0(72),
 336   fCalRocArrayT0Err(72),
 337   fCalRocArrayQ(72),
 338   fCalRocArrayRMS(72),
 339   fCalRocArrayOutliers(72),
 340   fHistoQArray(72),
 341   fHistoT0Array(72),
 342   fHistoRMSArray(72),
 343   fMeanT0rms(0),
 344   fMeanQrms(0),
 345   fMeanRMSrms(0),
 346   fHistoTmean(72),
 347   fParamArrayEventPol1(72),
 348   fParamArrayEventPol2(72),
 349   fTMeanArrayEvent(72),
 350   fQMeanArrayEvent(72),
 351   fVEventTime(1000),
 352   fVEventNumber(1000),
 353   fVTime0SideA(1000),
 354   fVTime0SideC(1000),
 355   fEventId(-1),
 356   fOldRunNumber(0),
 357   fPadTimesArrayEvent(72),
 358   fPadQArrayEvent(72),
 359   fPadRMSArrayEvent(72),
 360   fPadPedestalArrayEvent(72),
 361   fCurrentChannel(-1),
 362   fCurrentSector(-1),
 363   fCurrentRow(-1),
 364   fMaxPadSignal(-1),
 365   fMaxTimeBin(-1),
 366 //   fPadSignal(1024),
 367   fPadPedestal(0),
 368   fPadNoise(0),
 369   fVTime0Offset(72),
 370   fVTime0OffsetCounter(72),
 371   fVMeanQ(72),
 372   fVMeanQCounter(72),
 373   fCurrentCETimeRef(0),
 374   fProcessOld(kTRUE),
 375   fProcessNew(kFALSE),
 376   fAnalyseNew(kTRUE),
 377   fHnDrift(0x0),
 378   fArrHnDrift(100),
 379   fTimeBursts(100),
 380   fArrFitGraphs(0x0),
 381   fEventInBunch(0)
 382 {
 383   //
 384   // AliTPCSignal default constructor
 385   //
 386   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 387   fFirstTimeBin=650;
 388   fLastTimeBin=1030;
 389   fParam->Update();
 390   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 391   for (Int_t i=0;i<14;++i){
 392     fPeaks[i]=0;
 393     fPeakWidths[i]=0;
 394   }
 395   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
 396 }
 397 //_____________________________________________________________________
 398 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
 399   AliTPCCalibRawBase(sig),
 400   fNbinsT0(sig.fNbinsT0),
 401   fXminT0(sig.fXminT0),
 402   fXmaxT0(sig.fXmaxT0),
 403   fNbinsQ(sig.fNbinsQ),
 404   fXminQ(sig.fXminQ),
 405   fXmaxQ(sig.fXmaxQ),
 406   fNbinsRMS(sig.fNbinsRMS),
 407   fXminRMS(sig.fXminRMS),
 408   fXmaxRMS(sig.fXmaxRMS),
 409   fPeakDetMinus(sig.fPeakDetMinus),
 410   fPeakDetPlus(sig.fPeakDetPlus),
 411   fPeakIntMinus(sig.fPeakIntMinus),
 412   fPeakIntPlus(sig.fPeakIntPlus),
 413   fNoiseThresholdMax(sig.fNoiseThresholdMax),
 414   fNoiseThresholdSum(sig.fNoiseThresholdSum),
 415   fIsZeroSuppressed(sig.fIsZeroSuppressed),
 416   fLastSector(-1),
 417   fSecRejectRatio(.4),
 418   fParam(new AliTPCParam),
 419   fPedestalTPC(0x0),
 420   fPadNoiseTPC(0x0),
 421   fPedestalROC(0x0),
 422   fPadNoiseROC(0x0),
 423   fCalRocArrayT0(72),
 424   fCalRocArrayT0Err(72),
 425   fCalRocArrayQ(72),
 426   fCalRocArrayRMS(72),
 427   fCalRocArrayOutliers(72),
 428   fHistoQArray(72),
 429   fHistoT0Array(72),
 430   fHistoRMSArray(72),
 431   fMeanT0rms(sig.fMeanT0rms),
 432   fMeanQrms(sig.fMeanQrms),
 433   fMeanRMSrms(sig.fMeanRMSrms),
 434   fHistoTmean(72),
 435   fParamArrayEventPol1(72),
 436   fParamArrayEventPol2(72),
 437   fTMeanArrayEvent(72),
 438   fQMeanArrayEvent(72),
 439   fVEventTime(sig.fVEventTime),
 440   fVEventNumber(sig.fVEventNumber),
 441   fVTime0SideA(sig.fVTime0SideA),
 442   fVTime0SideC(sig.fVTime0SideC),
 443   fEventId(-1),
 444   fOldRunNumber(0),
 445   fPadTimesArrayEvent(72),
 446   fPadQArrayEvent(72),
 447   fPadRMSArrayEvent(72),
 448   fPadPedestalArrayEvent(72),
 449   fCurrentChannel(-1),
 450   fCurrentSector(-1),
 451   fCurrentRow(-1),
 452   fMaxPadSignal(-1),
 453   fMaxTimeBin(-1),
 454 //   fPadSignal(1024),
 455   fPadPedestal(0),
 456   fPadNoise(0),
 457   fVTime0Offset(72),
 458   fVTime0OffsetCounter(72),
 459   fVMeanQ(72),
 460   fVMeanQCounter(72),
 461   fCurrentCETimeRef(0),
 462   fProcessOld(sig.fProcessOld),
 463   fProcessNew(sig.fProcessNew),
 464   fAnalyseNew(sig.fAnalyseNew),
 465   fHnDrift(0x0),
 466   fArrHnDrift(100),
 467   fTimeBursts(100),
 468   fArrFitGraphs(0x0),
 469   fEventInBunch(0)
 470 {
 471   //
 472   // AliTPCSignal copy constructor
 473   //
 474   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 475
 476   for (Int_t iSec = 0; iSec < 72; ++iSec){
 477     const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 478     const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 479     const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 480     const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 481
 482     const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 483     const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 484     const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 485
 486     if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 487     if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 488     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 489     if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 490
 491     if ( hQ != 0x0 ){
 492       TH2S *hNew = new TH2S(*hQ);
 493       hNew->SetDirectory(0);
 494       fHistoQArray.AddAt(hNew,iSec);
 495     }
 496     if ( hT0 != 0x0 ){
 497       TH2S *hNew = new TH2S(*hT0);
 498       hNew->SetDirectory(0);
 499       fHistoT0Array.AddAt(hNew,iSec);
 500     }
 501     if ( hRMS != 0x0 ){
 502       TH2S *hNew = new TH2S(*hRMS);
 503       hNew->SetDirectory(0);
 504       fHistoRMSArray.AddAt(hNew,iSec);
 505     }
 506   }
 507
 508   //copy fit parameters event by event
 509   TObjArray *arr=0x0;
 510   for (Int_t iSec=0; iSec<72; ++iSec){
 511     arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
 512     if ( arr ){
 513       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 514       fParamArrayEventPol1.AddAt(arrEvents, iSec);
 515       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 516         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 517           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 518     }
 519
 520     arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
 521     if ( arr ){
 522       TObjArray *arrEvents = new TObjArray(arr->GetSize());
 523       fParamArrayEventPol2.AddAt(arrEvents, iSec);
 524       for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
 525         if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
 526           arrEvents->AddAt(new TVectorD(*vec),iEvent);
 527     }
 528
 529     TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
 530     TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
 531     if ( vMeanTime )
 532       fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
 533     if ( vMeanQ )
 534       fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
 535   }
 536
 537
 538   fVEventTime.ResizeTo(sig.fVEventTime);
 539   fVEventNumber.ResizeTo(sig.fVEventNumber);
 540   fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
 541   fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
 542
 543   fParam->Update();
 544
 545   for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
 546     TObject *o=sig.fArrHnDrift.UncheckedAt(i);
 547     if (o){
 548       TObject *newo=o->Clone("fHnDrift");
 549       fArrHnDrift.AddAt(newo,i);
 550       if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
 551     }
 552   }
 553
 554   for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
 555     fTimeBursts[i]=sig.fTimeBursts[i];
 556   }
 557
 558   for (Int_t i=0;i<14;++i){
 559     fPeaks[i]=sig.fPeaks[i];
 560     fPeakWidths[i]=sig.fPeakWidths[i];
 561   }
 562   if (sig.fArrFitGraphs) {
 563     fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
 564     fArrFitGraphs->SetOwner();
 565   }
 566
 567   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
 568
 569 }
 570 //_____________________________________________________________________
 571 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
 572   AliTPCCalibRawBase(),
 573   fNbinsT0(200),
 574   fXminT0(-5),
 575   fXmaxT0(5),
 576   fNbinsQ(200),
 577   fXminQ(1),
 578   fXmaxQ(40),
 579   fNbinsRMS(100),
 580   fXminRMS(0.1),
 581   fXmaxRMS(5.1),
 582   fPeakDetMinus(2),
 583   fPeakDetPlus(3),
 584   fPeakIntMinus(2),
 585   fPeakIntPlus(2),
 586   fNoiseThresholdMax(5.),
 587   fNoiseThresholdSum(8.),
 588   fIsZeroSuppressed(kFALSE),
 589   fLastSector(-1),
 590   fSecRejectRatio(.4),
 591   fParam(new  AliTPCParam),
 592   fPedestalTPC(0x0),
 593   fPadNoiseTPC(0x0),
 594   fPedestalROC(0x0),
 595   fPadNoiseROC(0x0),
 596   fCalRocArrayT0(72),
 597   fCalRocArrayT0Err(72),
 598   fCalRocArrayQ(72),
 599   fCalRocArrayRMS(72),
 600   fCalRocArrayOutliers(72),
 601   fHistoQArray(72),
 602   fHistoT0Array(72),
 603   fHistoRMSArray(72),
 604   fMeanT0rms(0),
 605   fMeanQrms(0),
 606   fMeanRMSrms(0),
 607   fHistoTmean(72),
 608   fParamArrayEventPol1(72),
 609   fParamArrayEventPol2(72),
 610   fTMeanArrayEvent(72),
 611   fQMeanArrayEvent(72),
 612   fVEventTime(1000),
 613   fVEventNumber(1000),
 614   fVTime0SideA(1000),
 615   fVTime0SideC(1000),
 616   fEventId(-1),
 617   fOldRunNumber(0),
 618   fPadTimesArrayEvent(72),
 619   fPadQArrayEvent(72),
 620   fPadRMSArrayEvent(72),
 621   fPadPedestalArrayEvent(72),
 622   fCurrentChannel(-1),
 623   fCurrentSector(-1),
 624   fCurrentRow(-1),
 625   fMaxPadSignal(-1),
 626   fMaxTimeBin(-1),
 627 //   fPadSignal(1024),
 628   fPadPedestal(0),
 629   fPadNoise(0),
 630   fVTime0Offset(72),
 631   fVTime0OffsetCounter(72),
 632   fVMeanQ(72),
 633   fVMeanQCounter(72),
 634   fCurrentCETimeRef(0),
 635   fProcessOld(kTRUE),
 636   fProcessNew(kFALSE),
 637   fAnalyseNew(kTRUE),
 638   fHnDrift(0x0),
 639   fArrHnDrift(100),
 640   fTimeBursts(100),
 641   fArrFitGraphs(0x0),
 642   fEventInBunch(0)
 643 {
 644   //
 645   // constructor which uses a tmap as input to set some specific parameters
 646   //
 647   SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
 648   fFirstTimeBin=650;
 649   fLastTimeBin=1030;
 650   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 651   if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 652   if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
 653   if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
 654   if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
 655   if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
 656   if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
 657   if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
 658   if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
 659   if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
 660   if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
 661   if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
 662   if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
 663   if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
 664   if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
 665   if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
 666   if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
 667   if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
 668   if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
 669   if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
 670
 671   if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
 672   if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
 673   if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
 674
 675   for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
 676   for (Int_t i=0;i<14;++i){
 677     fPeaks[i]=0;
 678     fPeakWidths[i]=0;
 679   }
 680
 681   fParam->Update();
 682   for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
 683 }
 684
 685 //_____________________________________________________________________
 686 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
 687 {
 688   //
 689   // assignment operator
 690   //
 691   if (&source == this) return *this;
 692   new (this) AliTPCCalibCE(source);
 693
 694   return *this;
 695 }
 696 //_____________________________________________________________________
 697 AliTPCCalibCE::~AliTPCCalibCE()
 698 {
 699   //
 700   // destructor
 701   //
 702
 703   fCalRocArrayT0.Delete();
 704   fCalRocArrayT0Err.Delete();
 705   fCalRocArrayQ.Delete();
 706   fCalRocArrayRMS.Delete();
 707   fCalRocArrayOutliers.Delete();
 708
 709   fHistoQArray.Delete();
 710   fHistoT0Array.Delete();
 711   fHistoRMSArray.Delete();
 712
 713   fHistoTmean.Delete();
 714
 715   fParamArrayEventPol1.Delete();
 716   fParamArrayEventPol2.Delete();
 717   fTMeanArrayEvent.Delete();
 718   fQMeanArrayEvent.Delete();
 719
 720   fPadTimesArrayEvent.Delete();
 721   fPadQArrayEvent.Delete();
 722   fPadRMSArrayEvent.Delete();
 723   fPadPedestalArrayEvent.Delete();
 724
 725   fArrHnDrift.SetOwner();
 726   fArrHnDrift.Delete();
 727
 728   if (fArrFitGraphs){
 729     fArrFitGraphs->SetOwner();
 730     delete fArrFitGraphs;
 731   }
 732 }
 733 //_____________________________________________________________________
 734 Int_t AliTPCCalibCE::Update(const Int_t icsector,
 735                                 const Int_t icRow,
 736                                 const Int_t icPad,
 737                                 const Int_t icTimeBin,
 738                                 const Float_t csignal)
 739 {
 740   //
 741   // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 742   // no extra analysis necessary. Assumes knowledge of the signal shape!
 743   // assumes that it is looped over consecutive time bins of one pad
 744   //
 745
 746   if (!fProcessOld) return 0;
 747   //temp
 748
 749   if (icRow<0) return 0;
 750   if (icPad<0) return 0;
 751   if (icTimeBin<0) return 0;
 752   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 753
 754   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 755
 756   //init first pad and sector in this event
 757   if ( fCurrentChannel == -1 ) {
 758     fLastSector=-1;
 759     fCurrentChannel = iChannel;
 760     fCurrentSector  = icsector;
 761     fCurrentRow     = icRow;
 762   }
 763
 764   //process last pad if we change to a new one
 765   if ( iChannel != fCurrentChannel ){
 766     ProcessPad();
 767     fLastSector=fCurrentSector;
 768     fCurrentChannel = iChannel;
 769     fCurrentSector  = icsector;
 770     fCurrentRow     = icRow;
 771   }
 772
 773   //fill signals for current pad
 774   fPadSignal[icTimeBin]=csignal;
 775   if ( csignal > fMaxPadSignal ){
 776     fMaxPadSignal = csignal;
 777     fMaxTimeBin   = icTimeBin;
 778   }
 779   return 0;
 780 }
 781
 782 //_____________________________________________________________________
 783 void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
 784                   const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
 785 {
 786   //
 787   // new filling method to fill the THnSparse histogram
 788   //
 789
 790   //only in new processing mode
 791   if (!fProcessNew) return;
 792     //don't use the IROCs and inner part of OROCs
 793   if (sector<36) return;
 794   if (row<40) return;
 795   //only bunches with reasonable length
 796   if (length<3||length>10) return;
 797
 798   UShort_t  timeBin    = (UShort_t)startTimeBin;
 799   //skip first laser layer
 800   if (timeBin<280) return;
 801
 802   Double_t timeBurst=SetBurstHnDrift();
 803
 804   Int_t cePeak=((sector/18)%2)*7+6;
 805   //after 1 event setup peak ranges
 806   if (fEventInBunch==1 && fPeaks[cePeak]==0) {
 807     // set time range
 808     fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
 809     FindLaserLayers();
 810     // set time range
 811     fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
 812     fHnDrift->Reset();
 813   }
 814
 815   // After the first event only fill every 5th  bin in a row with the CE information
 816   Int_t padFill=pad;
 817   if (fEventInBunch==0||(fPeaks[cePeak]>100&&TMath::Abs((Short_t)fPeaks[cePeak]-(Short_t)timeBin)<(Short_t)fPeakWidths[cePeak])){
 818     Int_t mod=5;
 819     Int_t n=pad/mod;
 820     padFill=mod*n+mod/2;
 821   }
 822
 823   //noise removal
 824   if (!IsPeakInRange(timeBin+length/2,sector)) return;
 825
 826   Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
 827       (Double_t)padFill,(Double_t)timeBin,timeBurst};
 828
 829   for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
 830     Float_t sig=(Float_t)signal[iTimeBin];
 831      // if (fPeaks[6]>900&&timeBin>(fPeaks[6]-20)&&sig<20) continue;
 832      // if (fPeaks[6]>900&&timeBin<(fPeaks[6]-fPeakWidth[6])&&sig<5) continue;
 833     x[3]=timeBin;
 834     fHnDrift->Fill(x,sig);
 835     --timeBin;
 836   }
 837 }
 838 //_____________________________________________________________________
 839 void AliTPCCalibCE::FindLaserLayers()
 840 {
 841   //
 842   // Find the laser layer positoins
 843   //
 844
 845   //A-side + C-side
 846   for (Int_t iside=0;iside<2;++iside){
 847     Int_t add=7*iside;
 848     //find CE signal position and width
 849     fHnDrift->GetAxis(0)->SetRangeUser(36+iside*18,53+iside*18);
 850     TH1D *hproj=fHnDrift->Projection(3);
 851     hproj->GetXaxis()->SetRangeUser(700,1030);
 852     Int_t maxbin=hproj->GetMaximumBin();
 853     Double_t binc=hproj->GetBinCenter(maxbin);
 854     hproj->GetXaxis()->SetRangeUser(binc-5,binc+5);
 855
 856     fPeaks[add+6]=(UShort_t)TMath::Nint(binc);
 857   //   fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
 858     fPeakWidths[add+6]=7;
 859
 860     hproj->GetXaxis()->SetRangeUser(0,maxbin-10);
 861     TSpectrum s(6);
 862     s.Search(hproj,2,"goff");
 863     Int_t index[6];
 864     TMath::Sort(6,s.GetPositionX(),index,kFALSE);
 865     for (Int_t i=0; i<6; ++i){
 866       fPeaks[i+add]=(UShort_t)TMath::Nint(s.GetPositionX()[index[i]]);
 867       fPeakWidths[i+add]=5;
 868     }
 869
 870     //other peaks
 871
 872 //    Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
 873 //    Int_t width=100;
 874
 875 //    for (Int_t i=3; i>=0; --i){
 876 //      hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
 877 //      fPeaks[i]=hproj->GetMaximumBin();
 878 //      fPeakWidths[i]=(UShort_t)TMath::Nint(10.);
 879 //      width=250;
 880 //      timepos=fPeaks[i]-width/2;
 881 //    }
 882
 883 //    for (Int_t i=add; i<add+7; ++i){
 884 //      printf("Peak: %u +- %u\n",fPeaks[i],fPeakWidths[i]);
 885 //    }
 886     //check width and reset peak if >100
 887 //    for (Int_t i=0; i<5; ++i){
 888 //      if (fPeakWidths[i]>100) {
 889 //        fPeaks[i]=0;
 890 //        fPeakWidths[i]=0;
 891 //      }
 892 //    }
 893
 894     delete hproj;
 895   }
 896 }
 897
 898 //_____________________________________________________________________
 899 void AliTPCCalibCE::FindPedestal(Float_t part)
 900 {
 901   //
 902     // find pedestal and noise for the current pad. Use either database or
 903     // truncated mean with part*100%
 904   //
 905   Bool_t noPedestal = kTRUE;
 906
 907     //use pedestal database if set
 908   if (fPedestalTPC&&fPadNoiseTPC){
 909         //only load new pedestals if the sector has changed
 910     if ( fCurrentSector!=fLastSector ){
 911       fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 912       fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 913     }
 914
 915     if ( fPedestalROC&&fPadNoiseROC ){
 916       fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
 917       fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 918       noPedestal   = kFALSE;
 919     }
 920
 921   }
 922
 923     //if we are not running with pedestal database, or for the current sector there is no information
 924     //available, calculate the pedestal and noise on the fly
 925   if ( noPedestal ) {
 926     fPadPedestal = 0;
 927     fPadNoise    = 0;
 928     if ( fIsZeroSuppressed ) return;
 929     const Int_t kPedMax = 100;  //maximum pedestal value
 930     Float_t  max    =  0;
 931     Float_t  maxPos =  0;
 932     Int_t    median =  -1;
 933     Int_t    count0 =  0;
 934     Int_t    count1 =  0;
 935     //
 936     Float_t padSignal=0;
 937     //
 938     UShort_t histo[kPedMax];
 939     memset(histo,0,kPedMax*sizeof(UShort_t));
 940
 941         //fill pedestal histogram
 942     for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 943       padSignal = fPadSignal[i];
 944       if (padSignal<=0) continue;
 945       if (padSignal>max && i>10) {
 946         max = padSignal;
 947         maxPos = i;
 948       }
 949       if (padSignal>kPedMax-1) continue;
 950       histo[int(padSignal+0.5)]++;
 951       count0++;
 952     }
 953         //find median
 954     for (Int_t i=1; i<kPedMax; ++i){
 955       if (count1<count0*0.5) median=i;
 956       count1+=histo[i];
 957     }
 958         // truncated mean
 959     //
 960     Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 961     //
 962     for (Int_t idelta=1; idelta<10; ++idelta){
 963       if (median-idelta<=0) continue;
 964       if (median+idelta>kPedMax) continue;
 965       if (count<part*count1){
 966         count+=histo[median-idelta];
 967         mean +=histo[median-idelta]*(median-idelta);
 968         rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 969         count+=histo[median+idelta];
 970         mean +=histo[median+idelta]*(median+idelta);
 971         rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 972       }
 973     }
 974     if ( count > 0 ) {
 975       mean/=count;
 976       rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 977       fPadPedestal = mean;
 978       fPadNoise    = rms;
 979     }
 980   }
 981 }
 982 //_____________________________________________________________________
 983 void AliTPCCalibCE::UpdateCETimeRef()
 984 {
 985   // Find the time reference of the last valid CE signal in sector
 986   // for irocs of the A-Side the reference of the corresponging OROC is returned
 987   // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
 988   if ( fLastSector == fCurrentSector ) return;
 989   Int_t sector=fCurrentSector;
 990   if ( sector < 18 ) sector+=36;
 991   fCurrentCETimeRef=0;
 992   TVectorF *vtRef = GetTMeanEvents(sector);
 993   if ( !vtRef ) return;
 994   Int_t vtRefSize= vtRef->GetNrows();
 995   if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
 996   else vtRefSize=fNevents;
 997   while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
 998   fCurrentCETimeRef=(*vtRef)[vtRefSize];
 999   AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
1000 }
1001 //_____________________________________________________________________
1002 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
1003 {
1004   //
1005     //  Find position, signal width and height of the CE signal (last signal)
1006     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
1007     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
1008   //
1009
1010   Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
1011   Int_t   cemaxpos       = 0;
1012   Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;  // threshold for the signal sum
1013   const Int_t    kCemin  = fPeakIntMinus;             // range for the analysis of the ce signal +- channels from the peak
1014   const Int_t    kCemax  = fPeakIntPlus;
1015
1016   Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
1017
1018     // find maximum closest to the sector mean from the last event
1019   for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
1020         // get sector mean of last event
1021     Float_t tmean = fCurrentCETimeRef;
1022     if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
1023       minDist  = tmean-maxima[imax];
1024       cemaxpos = (Int_t)maxima[imax];
1025     }
1026   }
1027 //   printf("L1 phase TB: %f\n",GetL1PhaseTB());
1028   if (cemaxpos!=0){
1029     ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
1030     for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
1031       if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
1032         Float_t signal = fPadSignal[i]-fPadPedestal;
1033         if (signal>0) {
1034           ceTime+=signal*(i+0.5);
1035           ceRMS +=signal*(i+0.5)*(i+0.5);
1036           ceQsum+=signal;
1037         }
1038       }
1039     }
1040   }
1041   if (ceQmax&&ceQsum>ceSumThreshold) {
1042     ceTime/=ceQsum;
1043     ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
1044     ceTime-=GetL1PhaseTB();
1045     fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
1046     fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1047
1048   //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1049   //                                the pick-up signal should scale with the pad area. In addition
1050   //                                the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1051   //                                ratio 2/3. The pad area we express in cm2. We normalise the signal
1052   //                                to the OROC signal (factor 2/3 for the IROCs).
1053     Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1054     if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1055
1056     ceQsum/=norm;
1057     fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1058     fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1059   } else {
1060     ceQmax=0;
1061     ceTime=0;
1062     ceRMS =0;
1063     ceQsum=0;
1064   }
1065   param[0] = ceQmax;
1066   param[1] = ceTime;
1067   param[2] = ceRMS;
1068   qSum     = ceQsum;
1069 }
1070 //_____________________________________________________________________
1071 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
1072 {
1073   //
1074   // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1075   // and 'tplus' timebins after 'pos'
1076   //
1077   if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1078   for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1079     if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1080   for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1081     if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1082       ++iTime2;
1083     if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1084   }
1085   return kTRUE;
1086 }
1087 //_____________________________________________________________________
1088 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1089 {
1090   //
1091     // Find local maxima on the pad signal and Histogram them
1092   //
1093   Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal
1094   Int_t   count       = 0;
1095
1096   for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1097     if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1098     if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
1099       if (count<maxima.GetNrows()){
1100         maxima.GetMatrixArray()[count++]=i;
1101         GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
1102         i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin  fPeakDetMinus+fPeakDetPlus-1
1103       }
1104     }
1105   }
1106 }
1107 //_____________________________________________________________________
1108 void AliTPCCalibCE::ProcessPad()
1109 {
1110   //
1111   //  Process data of current pad
1112   //
1113   FindPedestal();
1114
1115   TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
1116                              // + central electrode and possibly post peaks from the CE signal
1117                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
1118   FindLocalMaxima(maxima);
1119   if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
1120
1121   UpdateCETimeRef();                       // update the time refenrence for the current sector
1122   if ( fCurrentCETimeRef<1e-30 ) return;      //return if we don't have time 0 info, eg if only one side has laser
1123   TVectorD param(3);
1124   Float_t  qSum;
1125   FindCESignal(param, qSum, maxima);
1126
1127   Double_t meanT  = param[1];
1128   Double_t sigmaT = param[2];
1129
1130     //Fill Event T0 counter
1131   (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1132
1133     //Fill Q histogram
1134   GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1135
1136     //Fill RMS histogram
1137   GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1138
1139
1140     //Fill debugging info
1141   if ( GetStreamLevel()>0 ){
1142     (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
1143     (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
1144     (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
1145   }
1146
1147   ResetPad();
1148 }
1149 //_____________________________________________________________________
1150 void AliTPCCalibCE::EndEvent()
1151 {
1152   //  Process data of current pad
1153   //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
1154   //  before the EndEvent function to set the event timestamp and number!!!
1155   //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
1156   //  function was called
1157   if (!fProcessOld) {
1158     if (fProcessNew){
1159       ++fNevents;
1160       ++fEventInBunch;
1161     }
1162     return;
1163   }
1164
1165   //check if last pad has allready been processed, if not do so
1166   if ( fMaxTimeBin>-1 ) ProcessPad();
1167
1168   AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
1169
1170   TVectorD param(3);
1171   TMatrixD dummy(3,3);
1172 //    TVectorF vMeanTime(72);
1173 //    TVectorF vMeanQ(72);
1174   AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1175   AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1176
1177   //find mean time0 offset for side A and C
1178   //use only orocs due to the better statistics
1179   Double_t time0Side[2];       //time0 for side A:0 and C:1
1180   Double_t time0SideCount[2];  //time0 counter for side A:0 and C:1
1181   time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1182   for ( Int_t iSec = 36; iSec<72; ++iSec ){
1183     time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1184     time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1185   }
1186   if ( time0SideCount[0] >0  )
1187     time0Side[0]/=time0SideCount[0];
1188   if ( time0SideCount[1] >0 )
1189     time0Side[1]/=time0SideCount[1];
1190     // end find time0 offset
1191   AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1192   Int_t nSecMeanT=0;
1193   //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1194   for ( Int_t iSec = 0; iSec<72; ++iSec ){
1195     AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1196     //find median and then calculate the mean around it
1197     TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
1198     if ( !hMeanT ) continue;
1199     //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1200     if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1201       hMeanT->Reset();
1202       AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1203       continue;
1204     }
1205
1206     Double_t entries = hMeanT->GetEffectiveEntries();
1207     Double_t sum     = 0;
1208     Short_t *arr     = hMeanT->GetArray()+1;
1209     Int_t ibin=0;
1210     for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1211       sum+=arr[ibin];
1212       if ( sum>=(entries/2.) ) break;
1213     }
1214     Int_t delta = 4;
1215     Int_t firstBin = fFirstTimeBin+ibin-delta;
1216     Int_t lastBin  = fFirstTimeBin+ibin+delta;
1217     if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1218     if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
1219     Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1220
1221         // check boundaries for ebye info of mean time
1222     TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1223     Int_t vSize=vMeanTime->GetNrows();
1224     if ( vSize < fNevents+1 ){
1225       vMeanTime->ResizeTo(vSize+100);
1226     }
1227
1228     // store mean time for the readout sides
1229     vSize=fVTime0SideA.GetNrows();
1230     if ( vSize < fNevents+1 ){
1231       fVTime0SideA.ResizeTo(vSize+100);
1232       fVTime0SideC.ResizeTo(vSize+100);
1233     }
1234     fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1235     fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1236
1237     vMeanTime->GetMatrixArray()[fNevents]=median;
1238     nSecMeanT++;
1239     // end find median
1240
1241     TVectorF *vTimes = GetPadTimesEvent(iSec);
1242     if ( !vTimes ) continue;                     //continue if no time information for this sector is available
1243
1244     AliTPCCalROC calIrocOutliers(0);
1245     AliTPCCalROC calOrocOutliers(36);
1246
1247     // calculate mean Q of the sector
1248     TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1249     vSize=vMeanQ->GetNrows();
1250     if ( vSize < fNevents+1 ){
1251       vMeanQ->ResizeTo(vSize+100);
1252     }
1253     Float_t meanQ = 0;
1254     if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1255     vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1256
1257     for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1258       Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
1259
1260             //set values for temporary roc calibration class
1261       if ( iSec < 36 ) {
1262         calIroc->SetValue(iChannel, time);
1263         if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1264
1265       } else {
1266         calOroc->SetValue(iChannel, time);
1267         if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1268       }
1269
1270       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1271         // test that we really found the CE signal reliably
1272         if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1273           GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1274
1275
1276
1277             //-------------------------------  Debug start  ------------------------------
1278       if ( GetStreamLevel()>0 ){
1279         TTreeSRedirector *streamer=GetDebugStreamer();
1280         if (streamer){
1281           Int_t row=0;
1282           Int_t pad=0;
1283           Int_t padc=0;
1284
1285           Float_t q   = (*GetPadQEvent(iSec))[iChannel];
1286           Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1287
1288           UInt_t channel=iChannel;
1289           Int_t sector=iSec;
1290
1291           while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1292           pad = channel-fROC->GetRowIndexes(sector)[row];
1293           padc = pad-(fROC->GetNPads(sector,row)/2);
1294
1295 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1296 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
1297 //                                  fLastTimeBin-fFirstTimeBin,
1298 //                                  fFirstTimeBin,fLastTimeBin);
1299 //              h1->SetDirectory(0);
1300         //
1301 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1302 //                  h1->Fill(i,fPadSignal(i));
1303
1304           Double_t t0Sec = 0;
1305           if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1306             t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1307           Double_t t0Side = time0Side[(iSec/18)%2];
1308           (*streamer) << "DataPad" <<
1309             "Event=" << fNevents <<
1310             "RunNumber=" << fRunNumber <<
1311             "TimeStamp="   << fTimeStamp <<
1312             "Sector="<< sector <<
1313             "Row="   << row<<
1314             "Pad="   << pad <<
1315             "PadC="  << padc <<
1316             "PadSec="<< channel <<
1317             "Time0Sec="  << t0Sec <<
1318             "Time0Side=" << t0Side <<
1319             "Time="  << time <<
1320             "RMS="   << rms <<
1321             "Sum="   << q <<
1322             "MeanQ=" << meanQ <<
1323         //                  "hist.=" << h1 <<
1324             "\n";
1325
1326     //          delete h1;
1327         }
1328       }
1329       //-----------------------------  Debug end  ------------------------------
1330     }// end channel loop
1331
1332
1333     //do fitting now only in debug mode
1334     if (GetDebugLevel()>0){
1335       TVectorD paramPol1(3);
1336       TVectorD paramPol2(6);
1337       TMatrixD matPol1(3,3);
1338       TMatrixD matPol2(6,6);
1339       Float_t  chi2Pol1=0;
1340       Float_t  chi2Pol2=0;
1341
1342       if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1343         if ( iSec < 36 ){
1344           calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1345           calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1346         } else {
1347           calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1348           calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1349         }
1350
1351         GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1352         GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1353       }
1354
1355   //-------------------------------  Debug start  ------------------------------
1356       if ( GetStreamLevel()>0 ){
1357         TTreeSRedirector *streamer=GetDebugStreamer();
1358         if ( streamer ) {
1359           (*streamer) << "DataRoc" <<
1360 //    "Event=" << fEvent <<
1361             "RunNumber=" << fRunNumber <<
1362             "TimeStamp="   << fTimeStamp <<
1363             "Sector="<< iSec <<
1364             "hMeanT.=" << hMeanT <<
1365             "median=" << median <<
1366             "paramPol1.=" << &paramPol1 <<
1367             "paramPol2.=" << &paramPol2 <<
1368             "matPol1.="   << &matPol1 <<
1369             "matPol2.="   << &matPol2 <<
1370             "chi2Pol1="   << chi2Pol1 <<
1371             "chi2Pol2="   << chi2Pol2 <<
1372             "\n";
1373         }
1374       }
1375     }
1376         //-------------------------------  Debug end  ------------------------------
1377     hMeanT->Reset();
1378   }// end sector loop
1379     //return if no sector has a valid mean time
1380   if ( nSecMeanT == 0 ) return;
1381
1382
1383 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1384 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1385   if ( fVEventTime.GetNrows() < fNevents+1 ) {
1386     fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1387     fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1388   }
1389   fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1390   fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1391
1392   ++fNevents;
1393   if (fProcessNew) ++fEventInBunch;
1394   fOldRunNumber = fRunNumber;
1395
1396   delete calIroc;
1397   delete calOroc;
1398   AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1399 }
1400 //_____________________________________________________________________
1401 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1402                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1403                                   const Char_t *type, Bool_t force)
1404 {
1405     //
1406     // return pointer to TH2S histogram of 'type'
1407     // if force is true create a new histogram if it doesn't exist allready
1408     //
1409     if ( !force || arr->UncheckedAt(sector) )
1410       return (TH2S*)arr->UncheckedAt(sector);
1411
1412     // if we are forced and histogram doesn't exist yet create it
1413     // new histogram with Q calib information. One value for each pad!
1414     TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1415                           nbinsY, ymin, ymax,
1416                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1417     hist->SetDirectory(0);
1418     arr->AddAt(hist,sector);
1419     return hist;
1420 }
1421 //_____________________________________________________________________
1422 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1423 {
1424     //
1425     // return pointer to T0 histogram
1426     // if force is true create a new histogram if it doesn't exist allready
1427     //
1428     TObjArray *arr = &fHistoT0Array;
1429     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1430 }
1431 //_____________________________________________________________________
1432 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1433 {
1434     //
1435     // return pointer to Q histogram
1436     // if force is true create a new histogram if it doesn't exist allready
1437     //
1438     TObjArray *arr = &fHistoQArray;
1439     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1440 }
1441 //_____________________________________________________________________
1442 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1443 {
1444     //
1445     // return pointer to Q histogram
1446     // if force is true create a new histogram if it doesn't exist allready
1447     //
1448     TObjArray *arr = &fHistoRMSArray;
1449     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1450 }
1451 //_____________________________________________________________________
1452 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1453                               const Char_t *type, Bool_t force)
1454 {
1455     //
1456     // return pointer to TH1S histogram
1457     // if force is true create a new histogram if it doesn't exist allready
1458     //
1459     if ( !force || arr->UncheckedAt(sector) )
1460       return (TH1S*)arr->UncheckedAt(sector);
1461
1462     // if we are forced and histogram doesn't yes exist create it
1463     // new histogram with calib information. One value for each pad!
1464     TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1465                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1466     hist->SetDirectory(0);
1467     arr->AddAt(hist,sector);
1468     return hist;
1469 }
1470 //_____________________________________________________________________
1471 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1472 {
1473     //
1474     // return pointer to Q histogram
1475     // if force is true create a new histogram if it doesn't exist allready
1476     //
1477     TObjArray *arr = &fHistoTmean;
1478     return GetHisto(sector, arr, "LastTmean", force);
1479 }
1480 //_____________________________________________________________________
1481 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1482 {
1483     //
1484     // return pointer to Pad Info from 'arr' for the current event and sector
1485     // if force is true create it if it doesn't exist allready
1486     //
1487     if ( !force || arr->UncheckedAt(sector) )
1488         return (TVectorF*)arr->UncheckedAt(sector);
1489
1490     TVectorF *vect = new TVectorF(size);
1491     arr->AddAt(vect,sector);
1492     return vect;
1493 }
1494 //_____________________________________________________________________
1495 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1496 {
1497     //
1498     // return pointer to Pad Times Array for the current event and sector
1499     // if force is true create it if it doesn't exist allready
1500     //
1501     TObjArray *arr = &fPadTimesArrayEvent;
1502     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1503 }
1504 //_____________________________________________________________________
1505 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1506 {
1507     //
1508     // return pointer to Pad Q Array for the current event and sector
1509     // if force is true create it if it doesn't exist allready
1510     // for debugging purposes only
1511     //
1512
1513     TObjArray *arr = &fPadQArrayEvent;
1514     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1515 }
1516 //_____________________________________________________________________
1517 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1518 {
1519     //
1520     // return pointer to Pad RMS Array for the current event and sector
1521     // if force is true create it if it doesn't exist allready
1522     // for debugging purposes only
1523     //
1524     TObjArray *arr = &fPadRMSArrayEvent;
1525     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1526 }
1527 //_____________________________________________________________________
1528 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1529 {
1530     //
1531     // return pointer to Pad RMS Array for the current event and sector
1532     // if force is true create it if it doesn't exist allready
1533     // for debugging purposes only
1534     //
1535     TObjArray *arr = &fPadPedestalArrayEvent;
1536     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1537 }
1538 //_____________________________________________________________________
1539 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1540 {
1541     //
1542     // return pointer to the EbyE info of the mean arrival time for 'sector'
1543     // if force is true create it if it doesn't exist allready
1544     //
1545     TObjArray *arr = &fTMeanArrayEvent;
1546     return GetVectSector(sector,arr,100,force);
1547 }
1548 //_____________________________________________________________________
1549 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1550 {
1551     //
1552     // return pointer to the EbyE info of the mean arrival time for 'sector'
1553     // if force is true create it if it doesn't exist allready
1554     //
1555     TObjArray *arr = &fQMeanArrayEvent;
1556     return GetVectSector(sector,arr,100,force);
1557 }
1558 //_____________________________________________________________________
1559 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1560 {
1561     //
1562     // return pointer to ROC Calibration
1563     // if force is true create a new histogram if it doesn't exist allready
1564     //
1565     if ( !force || arr->UncheckedAt(sector) )
1566         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1567
1568     // if we are forced and histogram doesn't yes exist create it
1569
1570     // new AliTPCCalROC for T0 information. One value for each pad!
1571     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1572     arr->AddAt(croc,sector);
1573     return croc;
1574 }
1575 //_____________________________________________________________________
1576 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1577 {
1578     //
1579     // return pointer to Time 0 ROC Calibration
1580     // if force is true create a new histogram if it doesn't exist allready
1581     //
1582     TObjArray *arr = &fCalRocArrayT0;
1583     return GetCalRoc(sector, arr, force);
1584 }
1585 //_____________________________________________________________________
1586 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1587 {
1588     //
1589     // return pointer to the error of Time 0 ROC Calibration
1590     // if force is true create a new histogram if it doesn't exist allready
1591     //
1592     TObjArray *arr = &fCalRocArrayT0Err;
1593     return GetCalRoc(sector, arr, force);
1594 }
1595 //_____________________________________________________________________
1596 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1597 {
1598     //
1599     // return pointer to T0 ROC Calibration
1600     // if force is true create a new histogram if it doesn't exist allready
1601     //
1602     TObjArray *arr = &fCalRocArrayQ;
1603     return GetCalRoc(sector, arr, force);
1604 }
1605 //_____________________________________________________________________
1606 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1607 {
1608     //
1609     // return pointer to signal width ROC Calibration
1610     // if force is true create a new histogram if it doesn't exist allready
1611     //
1612     TObjArray *arr = &fCalRocArrayRMS;
1613     return GetCalRoc(sector, arr, force);
1614 }
1615 //_____________________________________________________________________
1616 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1617 {
1618     //
1619     // return pointer to Outliers
1620     // if force is true create a new histogram if it doesn't exist allready
1621     //
1622     TObjArray *arr = &fCalRocArrayOutliers;
1623     return GetCalRoc(sector, arr, force);
1624 }
1625 //_____________________________________________________________________
1626 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1627 {
1628     //
1629     // return pointer to TObjArray of fit parameters
1630     // if force is true create a new histogram if it doesn't exist allready
1631     //
1632     if ( !force || arr->UncheckedAt(sector) )
1633         return (TObjArray*)arr->UncheckedAt(sector);
1634
1635     // if we are forced and array doesn't yes exist create it
1636
1637     // new TObjArray for parameters
1638     TObjArray *newArr = new TObjArray;
1639     arr->AddAt(newArr,sector);
1640     return newArr;
1641 }
1642 //_____________________________________________________________________
1643 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1644 {
1645     //
1646     // return pointer to TObjArray of fit parameters from plane fit
1647     // if force is true create a new histogram if it doesn't exist allready
1648     //
1649     TObjArray *arr = &fParamArrayEventPol1;
1650     return GetParamArray(sector, arr, force);
1651 }
1652 //_____________________________________________________________________
1653 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1654 {
1655     //
1656     // return pointer to TObjArray of fit parameters from parabola fit
1657     // if force is true create a new histogram if it doesn't exist allready
1658     //
1659     TObjArray *arr = &fParamArrayEventPol2;
1660     return GetParamArray(sector, arr, force);
1661 }
1662
1663 //_____________________________________________________________________
1664 void AliTPCCalibCE::CreateDVhist()
1665 {
1666   //
1667   // Setup the THnSparse for the drift velocity determination
1668   //
1669
1670   //HnSparse bins
1671   //roc, row, pad, timebin, timestamp
1672   TTimeStamp begin(2010,01,01,0,0,0);
1673   TTimeStamp end(2035,01,01,0,0,0);
1674   Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1675
1676   Int_t    bins[kHnBinsDV] = { 72,  96,  140,  1030, nbinsTime};
1677   Double_t xmin[kHnBinsDV] = { 0.,  0.,   0.,    0., (Double_t)begin.GetSec()};
1678   Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1679
1680   fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1681   fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1682   fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1683   fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1684   fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1685   fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1686   fHnDrift->Reset();
1687 }
1688
1689 //_____________________________________________________________________
1690 void AliTPCCalibCE::ResetEvent()
1691 {
1692     //
1693     //  Reset global counters  -- Should be called before each event is processed
1694     //
1695     fLastSector=-1;
1696     fCurrentSector=-1;
1697     fCurrentRow=-1;
1698     fCurrentChannel=-1;
1699
1700     ResetPad();
1701
1702     fPadTimesArrayEvent.Delete();
1703     fPadQArrayEvent.Delete();
1704     fPadRMSArrayEvent.Delete();
1705     fPadPedestalArrayEvent.Delete();
1706
1707     for ( Int_t i=0; i<72; ++i ){
1708         fVTime0Offset.GetMatrixArray()[i]=0;
1709         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1710         fVMeanQ.GetMatrixArray()[i]=0;
1711         fVMeanQCounter.GetMatrixArray()[i]=0;
1712     }
1713 }
1714 //_____________________________________________________________________
1715 void AliTPCCalibCE::ResetPad()
1716 {
1717     //
1718     //  Reset pad infos -- Should be called after a pad has been processed
1719     //
1720     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1721       fPadSignal[i] = 0;
1722     fMaxTimeBin   = -1;
1723     fMaxPadSignal = -1;
1724     fPadPedestal  = -1;
1725     fPadNoise     = -1;
1726 }
1727 //_____________________________________________________________________
1728 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1729 {
1730   //
1731   //  Merge ce to the current AliTPCCalibCE
1732   //
1733   MergeBase(ce);
1734   Int_t nCEevents = ce->GetNeventsProcessed();
1735
1736   if (fProcessOld&&ce->fProcessOld){
1737   //merge histograms
1738     for (Int_t iSec=0; iSec<72; ++iSec){
1739       TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1740       TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1741       TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1742
1743
1744       if ( hRefQmerge ){
1745         TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1746         TH2S *hRefQ   = GetHistoQ(iSec);
1747         if ( hRefQ ) hRefQ->Add(hRefQmerge);
1748         else {
1749           TH2S *hist = new TH2S(*hRefQmerge);
1750           hist->SetDirectory(0);
1751           fHistoQArray.AddAt(hist, iSec);
1752         }
1753         hRefQmerge->SetDirectory(dir);
1754       }
1755       if ( hRefT0merge ){
1756         TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1757         TH2S *hRefT0  = GetHistoT0(iSec);
1758         if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1759         else {
1760           TH2S *hist = new TH2S(*hRefT0merge);
1761           hist->SetDirectory(0);
1762           fHistoT0Array.AddAt(hist, iSec);
1763         }
1764         hRefT0merge->SetDirectory(dir);
1765       }
1766       if ( hRefRMSmerge ){
1767         TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1768         TH2S *hRefRMS = GetHistoRMS(iSec);
1769         if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1770         else {
1771           TH2S *hist = new TH2S(*hRefRMSmerge);
1772           hist->SetDirectory(0);
1773           fHistoRMSArray.AddAt(hist, iSec);
1774         }
1775         hRefRMSmerge->SetDirectory(dir);
1776       }
1777
1778     }
1779
1780     // merge time information
1781
1782
1783     for (Int_t iSec=0; iSec<72; ++iSec){
1784       TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1785       TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1786       TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1787       TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1788
1789       TObjArray *arrPol1  = 0x0;
1790       TObjArray *arrPol2  = 0x0;
1791       TVectorF *vMeanTime = 0x0;
1792       TVectorF *vMeanQ    = 0x0;
1793
1794   //resize arrays
1795       if ( arrPol1CE && arrPol2CE ){
1796         arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1797         arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1798         arrPol1->Expand(fNevents+nCEevents);
1799         arrPol2->Expand(fNevents+nCEevents);
1800       }
1801       if ( vMeanTimeCE && vMeanQCE ){
1802         vMeanTime = GetTMeanEvents(iSec,kTRUE);
1803         vMeanQ    = GetQMeanEvents(iSec,kTRUE);
1804         vMeanTime->ResizeTo(fNevents+nCEevents);
1805         vMeanQ->ResizeTo(fNevents+nCEevents);
1806       }
1807
1808       for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1809         if ( arrPol1CE && arrPol2CE ){
1810           TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1811           TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1812           if ( paramPol1 && paramPol2 ){
1813             GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1814             GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1815           }
1816         }
1817         if ( vMeanTimeCE && vMeanQCE ){
1818           vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1819           vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1820         }
1821       }
1822     }
1823
1824
1825
1826     const TVectorD&  eventTimes  = ce->fVEventTime;
1827     const TVectorD&  eventIds    = ce->fVEventNumber;
1828     const TVectorF&  time0SideA  = ce->fVTime0SideA;
1829     const TVectorF&  time0SideC  = ce->fVTime0SideC;
1830     fVEventTime.ResizeTo(fNevents+nCEevents);
1831     fVEventNumber.ResizeTo(fNevents+nCEevents);
1832     fVTime0SideA.ResizeTo(fNevents+nCEevents);
1833     fVTime0SideC.ResizeTo(fNevents+nCEevents);
1834
1835     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1836       Double_t evTime     = eventTimes.GetMatrixArray()[iEvent];
1837       Double_t evId       = eventIds.GetMatrixArray()[iEvent];
1838       Float_t  t0SideA    = time0SideA.GetMatrixArray()[iEvent];
1839       Float_t  t0SideC    = time0SideC.GetMatrixArray()[iEvent];
1840
1841       fVEventTime.GetMatrixArray()[fNevents+iEvent]   = evTime;
1842       fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1843       fVTime0SideA.GetMatrixArray()[fNevents+iEvent]  = t0SideA;
1844       fVTime0SideC.GetMatrixArray()[fNevents+iEvent]  = t0SideC;
1845     }
1846   }
1847
1848   if (fProcessNew&&ce->fProcessNew) {
1849     if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1850       AliError("Number of bursts in the instances to merge are different. No merging done!");
1851     } else {
1852       for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1853         THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1854         THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1855         if (h && hce) h->Add(hce);
1856         else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1857       }
1858       //TODO: What to do with fTimeBursts???
1859     }
1860   }
1861
1862   fNevents+=nCEevents; //increase event counter
1863 }
1864
1865 //_____________________________________________________________________
1866 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1867 {
1868   //
1869   // Merge all objects of this type in list
1870   //
1871
1872   Long64_t nmerged=1;
1873
1874   TIter next(list);
1875   AliTPCCalibCE *ce=0;
1876   TObject *o=0;
1877
1878   while ( (o=next()) ){
1879     ce=dynamic_cast<AliTPCCalibCE*>(o);
1880     if (ce){
1881       Merge(ce);
1882       ++nmerged;
1883     }
1884   }
1885
1886   return nmerged;
1887 }
1888
1889 //_____________________________________________________________________
1890 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1891 {
1892   //
1893   // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1894   // or side (-1: A-Side, -2: C-Side)
1895   // xVariable:    0-event time, 1-event id, 2-internal event counter
1896   // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1897   // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1898   //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1899   //                               not used for mean time and mean Q )
1900   // for an example see class description at the beginning
1901   //
1902
1903   TVectorD *xVar = 0x0;
1904   TObjArray *aType = 0x0;
1905   Int_t npoints=0;
1906
1907   // sanity checks
1908   if ( (sector<-2)   || (sector>71)   ) return 0x0;  //sector outside valid range
1909   if ( (xVariable<0) || (xVariable>2) ) return 0x0;  //invalid x-variable
1910   if ( (fitType<0)   || (fitType>3)   ) return 0x0;  //invalid fit type
1911   if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1912   if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1913   if ( sector<0 && fitType!=2) return 0x0;  //for side wise information only mean time is available
1914
1915   if (sector>=0){
1916     if ( fitType==0 ){
1917       if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1918       aType = &fParamArrayEventPol1;
1919       if ( aType->At(sector)==0x0 ) return 0x0;
1920     }
1921     else if ( fitType==1 ){
1922       if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1923       aType = &fParamArrayEventPol2;
1924       if ( aType->At(sector)==0x0 ) return 0x0;
1925     }
1926
1927   }
1928   if ( xVariable == 0 ) xVar = &fVEventTime;
1929   if ( xVariable == 1 ) xVar = &fVEventNumber;
1930   if ( xVariable == 2 ) {
1931     xVar = new TVectorD(fNevents);
1932     for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1933   }
1934
1935   Double_t *x = new Double_t[fNevents];
1936   Double_t *y = new Double_t[fNevents];
1937
1938   for (Int_t ievent =0; ievent<fNevents; ++ievent){
1939     if ( fitType<2 ){
1940       TObjArray *events = (TObjArray*)(aType->At(sector));
1941       if ( events->GetSize()<=ievent ) break;
1942       TVectorD *v = (TVectorD*)(events->At(ievent));
1943       if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1944     } else if (fitType == 2) {
1945       Double_t xValue=(*xVar)[ievent];
1946       Double_t yValue=0;
1947       if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1948       else if (sector==-1) yValue=fVTime0SideA(ievent);
1949       else if (sector==-2) yValue=fVTime0SideC(ievent);
1950       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1951     }else if (fitType == 3) {
1952       Double_t xValue=(*xVar)[ievent];
1953       Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1954       if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1955     }
1956   }
1957
1958   TGraph *gr = new TGraph(npoints);
1959     //sort xVariable increasing
1960   Int_t    *sortIndex = new Int_t[npoints];
1961   TMath::Sort(npoints,x,sortIndex, kFALSE);
1962   for (Int_t i=0;i<npoints;++i){
1963     gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1964   }
1965
1966
1967   if ( xVariable == 2 ) delete xVar;
1968   delete [] x;
1969   delete [] y;
1970   delete [] sortIndex;
1971   return gr;
1972 }
1973 //_____________________________________________________________________
1974 void AliTPCCalibCE::Analyse()
1975 {
1976   //
1977   //  Calculate calibration constants
1978   //
1979
1980   if (fProcessOld){
1981     TVectorD paramQ(3);
1982     TVectorD paramT0(3);
1983     TVectorD paramRMS(3);
1984     TMatrixD dummy(3,3);
1985
1986     Float_t channelCounter=0;
1987     fMeanT0rms=0;
1988     fMeanQrms=0;
1989     fMeanRMSrms=0;
1990
1991     for (Int_t iSec=0; iSec<72; ++iSec){
1992       TH2S *hT0 = GetHistoT0(iSec);
1993       if (!hT0 ) continue;
1994
1995       AliTPCCalROC *rocQ     = GetCalRocQ  (iSec,kTRUE);
1996       AliTPCCalROC *rocT0    = GetCalRocT0 (iSec,kTRUE);
1997       AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1998       AliTPCCalROC *rocRMS   = GetCalRocRMS(iSec,kTRUE);
1999       AliTPCCalROC *rocOut   = GetCalRocOutliers(iSec,kTRUE);
2000
2001       TH2S *hQ   = GetHistoQ(iSec);
2002       TH2S *hRMS = GetHistoRMS(iSec);
2003
2004       Short_t *arrayhQ   = hQ->GetArray();
2005       Short_t *arrayhT0  = hT0->GetArray();
2006       Short_t *arrayhRMS = hRMS->GetArray();
2007
2008       UInt_t nChannels = fROC->GetNChannels(iSec);
2009
2010   //debug
2011       Int_t row=0;
2012       Int_t pad=0;
2013       Int_t padc=0;
2014   //! debug
2015
2016       for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
2017
2018
2019         Float_t cogTime0 = -1000;
2020         Float_t cogQ     = -1000;
2021         Float_t cogRMS   = -1000;
2022         Float_t cogOut   = 0;
2023         Float_t rms      = 0;
2024         Float_t rmsT0    = 0;
2025
2026
2027         Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
2028         Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
2029         Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
2030
2031         cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2032         fMeanQrms+=rms;
2033         cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2034         fMeanT0rms+=rmsT0;
2035         cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2036         fMeanRMSrms+=rms;
2037         channelCounter++;
2038
2039       /*
2040              //outlier specifications
2041          if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2042          cogOut = 1;
2043           cogTime0 = 0;
2044           cogQ     = 0;
2045           cogRMS   = 0;
2046       }
2047       */
2048         rocQ->SetValue(iChannel, cogQ*cogQ);
2049         rocT0->SetValue(iChannel, cogTime0);
2050         rocT0Err->SetValue(iChannel, rmsT0);
2051         rocRMS->SetValue(iChannel, cogRMS);
2052         rocOut->SetValue(iChannel, cogOut);
2053
2054
2055       //debug
2056         if ( GetStreamLevel() > 0 ){
2057           TTreeSRedirector *streamer=GetDebugStreamer();
2058           if ( streamer ) {
2059
2060             while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2061             pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2062             padc = pad-(fROC->GetNPads(iSec,row)/2);
2063
2064             (*streamer) << "DataEnd" <<
2065               "Sector="  << iSec      <<
2066               "Pad="     << pad       <<
2067               "PadC="    << padc      <<
2068               "Row="     << row       <<
2069               "PadSec="  << iChannel   <<
2070               "Q="       << cogQ      <<
2071               "T0="      << cogTime0  <<
2072               "RMS="     << cogRMS    <<
2073               "\n";
2074           }
2075         }
2076       //! debug
2077
2078       }
2079
2080     }
2081     if ( channelCounter>0 ){
2082       fMeanT0rms/=channelCounter;
2083       fMeanQrms/=channelCounter;
2084       fMeanRMSrms/=channelCounter;
2085     }
2086     //   if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2087     //    delete fDebugStreamer;
2088     //    fDebugStreamer = 0x0;
2089     fVEventTime.ResizeTo(fNevents);
2090     fVEventNumber.ResizeTo(fNevents);
2091     fVTime0SideA.ResizeTo(fNevents);
2092     fVTime0SideC.ResizeTo(fNevents);
2093   }
2094
2095   if (fProcessNew&&fAnalyseNew){
2096     AnalyseTrack();
2097     for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2098       THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2099       h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2100     }
2101   }
2102 }
2103
2104
2105
2106
2107 //
2108 // New functions that also use the laser tracks
2109 //
2110
2111
2112
2113 //_____________________________________________________________________
2114 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2115 {
2116   //
2117   //Find the local maximums for the projections to each axis
2118   //
2119
2120   //find laser layer positoins
2121   fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2122   FindLaserLayers();
2123   THnSparse *hProj=fHnDrift;
2124   Double_t posCE[4]={0.,0.,0.,0.};
2125   Double_t widthCE[4]={0.,0.,0.,0.};
2126
2127 //   if(fPeaks[4]!=0){
2128   // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2129
2130   for (Int_t i=0; i<4; ++i){
2131     Int_t ce=(i/2>0)*7+6;
2132     hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2133     TH1 *h=fHnDrift->Projection(3);
2134     h->GetXaxis()->SetRangeUser(fPeaks[ce]-fPeakWidths[ce],fPeaks[ce]+fPeakWidths[ce]);
2135     Int_t nbinMax=h->GetMaximumBin();
2136     Double_t maximum=h->GetMaximum();
2137 //     Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2138 //     if (nbinMax<700||maximum<maxExpected) continue;
2139     Double_t xbinMax=h->GetBinCenter(nbinMax);
2140     TF1 fgaus("gaus","gaus",xbinMax-5,xbinMax+5);
2141     fgaus.SetParameters(maximum,xbinMax,2);
2142     fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2143     fgaus.SetParLimits(2,0.2,4.);
2144     h->Fit(&fgaus,"RQN");
2145 //     Double_t deltaX=4*fgaus.GetParameter(2);
2146 //     xbinMax=fgaus.GetParameter(1);
2147     delete h;
2148     posCE[i]=fgaus.GetParameter(1);
2149     widthCE[i]=4*fgaus.GetParameter(2);
2150     hProj->GetAxis(0)->SetRangeUser(0,72);
2151   }
2152 //   }
2153   //Current drift velocity
2154   Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2155 //   cout<<"vdrift="<<vdrift<<endl;
2156
2157   AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2158   //loop over all entries in the histogram
2159   Int_t coord[5];
2160   for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2161     //get entry position and content
2162     Double_t adc=hProj->GetBinContent(ichunk,coord);
2163
2164
2165     Int_t sector = coord[0]-1;
2166     Int_t row    = coord[1]-1;
2167     Int_t pad    = coord[2]-1;
2168     Int_t timeBin= coord[3]-1;
2169     Double_t time   = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2170     Int_t side   = (sector/18)%2;
2171 //     return;
2172 //     fPeaks[4]=(UInt_t)posCE[sector/18];
2173 //     fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2174
2175     //cuts
2176     if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2177     if (adc < 5 ) continue;
2178     if (IsEdgePad(sector,row,pad)) continue;
2179 //     if (!IsPeakInRange(timeBin)) continue;
2180 //     if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2181 //     if (isCE&&(sector!=0)) continue;
2182
2183     Int_t padmin=-2, padmax=2;
2184     Int_t timemin=-2, timemax=2;
2185     Int_t minsumperpad=2;
2186     //CE or laser tracks
2187     Bool_t isCE=kFALSE;
2188     if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2189       isCE=kTRUE;
2190       padmin=0;
2191       padmax=0;
2192       timemin=-3;
2193       timemax=7;
2194     }
2195
2196     //
2197     // Find local maximum and cogs
2198     //
2199     Bool_t isMaximum=kTRUE;
2200     Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2201     Double_t cogY=0, rmsY=0;
2202     Int_t npart=0;
2203
2204     // for position calculation use
2205     for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2206       Float_t lxyz[3];
2207       fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2208
2209       for(Int_t itime=timemin;itime<=timemax;++itime){
2210
2211         Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2212         Double_t val=hProj->GetBinContent(a);
2213         totalmass+=val;
2214
2215         tbcm +=(timeBin+itime)*val;
2216         padcm+=(pad+ipad)*val;
2217         cogY +=lxyz[1]*val;
2218
2219         rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2220         rmspad+=(pad+ipad)*(pad+ipad)*val;
2221         rmsY  +=lxyz[1]*lxyz[1]*val;
2222
2223         if (val>0) ++npart;
2224         if (val>adc) {
2225           isMaximum=kFALSE;
2226           break;
2227         }
2228       }
2229       if (!isMaximum)  break;
2230     }
2231
2232     if (!isMaximum||!npart)  continue;
2233     if (totalmass<npart*minsumperpad) continue;
2234     if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2235                                     //for CE we want only one pad by construction
2236
2237     tbcm/=totalmass;
2238     padcm/=totalmass;
2239     cogY/=totalmass;
2240
2241     rmstb/=totalmass;
2242     rmspad/=totalmass;
2243     rmsY/=totalmass;
2244
2245     rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2246     rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2247     rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2248
2249     Int_t cog=TMath::Nint(padcm);
2250
2251     // timebin --> z position
2252     Float_t zlength=fParam->GetZLength(side);
2253 //     Float_t timePos=tbcm+(1000-fPeaks[4])
2254     // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2255     Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2256
2257     // local to global transformation--> x and y positions
2258     Float_t padlxyz[3];
2259     fROC->GetPositionLocal(sector,row,pad,padlxyz);
2260
2261     Double_t gxyz[3]={padlxyz[0],cogY,gz};
2262     Double_t lxyz[3]={padlxyz[0],cogY,gz};
2263     Double_t igxyz[3]={0,0,0};
2264     AliTPCTransform t1;
2265     t1.RotatedGlobal2Global(sector,gxyz);
2266
2267     Double_t mindist=0;
2268     Int_t trackID=-1;
2269     Int_t trackID2=-1;
2270
2271     //find track id and index of the position in the track (row)
2272     Int_t index=0;
2273     if (!isCE){
2274       index=row+(sector>35)*fROC->GetNRows(0);
2275       trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2276     } else {
2277       trackID=336+((sector/18)%2);
2278       index= fROC->GetRowIndexes(sector)[row]+pad; //  global pad position in sector
2279       if (sector<36) {
2280         index+=(sector%18)*fROC->GetNChannels(sector);
2281       } else {
2282         index+=18*fROC->GetNChannels(0);
2283         index+=(sector%18)*fROC->GetNChannels(sector);
2284       }
2285       //TODO: find out about the multiple peaks in the CE
2286 //       mindist=TMath::Abs(fPeaks[4]-tbcm);
2287       mindist=1.;
2288     }
2289
2290     // fill track vectors
2291     if (trackID>0){
2292       AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2293       Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2294
2295
2296 //      travel time effect of light includes
2297
2298       Double_t raylength=ltr->GetRayLength();
2299       Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2300       ltr->GetXYZ(globmir);
2301       if(trackID<336){
2302         if(side==0){
2303           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2304                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2305                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2306         }
2307         else {
2308           gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2309                                        +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2310                                        +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2311         }
2312       }
2313
2314       if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2315         ltr->fVecSec->GetMatrixArray()[index]=sector;
2316         ltr->fVecP2->GetMatrixArray()[index]=0;
2317         ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2318         ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2319         ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2320         ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2321         ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2322         ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2323         ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2324 //         ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2325       }
2326       TObjArray *arr=AliTPCLaserTrack::GetTracks();
2327       ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2328       igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2329       igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2330       igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2331     }
2332
2333
2334     if (fStreamLevel>4){
2335       (*GetDebugStreamer()) << "clusters" <<
2336         "run="   << fRunNumber <<
2337         "timestamp=" << timestamp <<
2338         "burst="     << burst     <<
2339         "side="      << side      <<
2340         "sec="       << sector    <<
2341         "row="       << row       <<
2342         "pad="       << pad       <<
2343         "padCog="    << cog       <<
2344         "timebin="   << timeBin   <<
2345         "cogCE="     << posCE[sector/18] <<
2346         "withCE="    << widthCE[sector/18] <<
2347         "index="     << index     <<
2348
2349         "padcm="     << padcm     <<
2350         "rmspad="    << rmspad    <<
2351
2352         "cogtb="     << tbcm      <<
2353         "rmstb="     << rmstb     <<
2354
2355         "npad="      << npart     <<
2356
2357         "lx="        << padlxyz[0]<<
2358         "ly="        << cogY      <<
2359         "lypad="     << padlxyz[1]<<
2360         "rmsY="      << rmsY      <<
2361
2362         "gx="        << gxyz[0]   <<
2363         "gy="        << gxyz[1]   <<
2364         "gz="        << gxyz[2]   <<
2365
2366         "igx="        << igxyz[0] <<
2367         "igy="        << igxyz[1] <<
2368         "igz="        << igxyz[2] <<
2369
2370         "mind="      << mindist   <<
2371         "max="       << adc       <<
2372         "trackid="   << trackID   <<
2373         "trackid2="   << trackID2   <<
2374         "npart="     << npart     <<
2375         "\n";
2376     } // end stream levelmgz.fElements
2377
2378   }
2379
2380 }
2381
2382 //_____________________________________________________________________
2383 void AliTPCCalibCE::AnalyseTrack()
2384 {
2385   //
2386   //  Analyse the tracks
2387   //
2388
2389
2390   AliTPCLaserTrack::LoadTracks();
2391 //   AliTPCParam *param=0x0;
2392 //   //cdb run number
2393 //   AliCDBManager *man=AliCDBManager::Instance();
2394 //   if (man->GetDefaultStorage()){
2395 //     AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2396 //     if (entry){
2397 //       entry->SetOwner(kTRUE);
2398 //       param = (AliTPCParam*)(entry->GetObject()->Clone());
2399 //     }
2400 //   }
2401 //   if (param){
2402 //     if (fParam) delete fParam;
2403 //     fParam=param;
2404 //   } else {
2405 //     AliError("Could not get updated AliTPCParam from OCDB!!!");
2406 //   }
2407
2408   //Measured and ideal laser tracks
2409   TObjArray* arrMeasured = SetupMeasured();
2410   TObjArray* arrIdeal    = AliTPCLaserTrack::GetTracks();
2411   AddCEtoIdeal(arrIdeal);
2412
2413   //find bursts and loop over them
2414   for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2415     Double_t timestamp=fTimeBursts[iburst];
2416     AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2417     fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2418     if (!fHnDrift) continue;
2419     UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2420     if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2421     fBinsLastAna[iburst]=entries;
2422
2423     for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2424 //     if (iburst==0) FindLaserLayers();
2425
2426     //reset laser tracks
2427     ResetMeasured(arrMeasured);
2428
2429     //find clusters and associate to the tracks
2430     FindLocalMaxima(arrMeasured, timestamp, iburst);
2431
2432     //calculate drift velocity
2433     CalculateDV(arrIdeal,arrMeasured,iburst);
2434
2435     //Dump information to file if requested
2436     if (fStreamLevel>2){
2437       //printf("make tree\n");
2438       //laser track information
2439
2440       for (Int_t itrack=0; itrack<338; ++itrack){
2441         TObject *iltr=arrIdeal->UncheckedAt(itrack);
2442         TObject *mltr=arrMeasured->UncheckedAt(itrack);
2443         (*GetDebugStreamer()) << "tracks" <<
2444           "run="   << fRunNumber <<
2445           "time=" << timestamp <<
2446           "burst="<< iburst <<
2447           "iltr.=" << iltr <<
2448           "mltr.=" << mltr <<
2449           "\n";
2450       }
2451     }
2452   }
2453   if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2454 }
2455
2456 //_____________________________________________________________________
2457 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2458                                       const Double_t *peakposloc, Int_t &itrackMin2)
2459 {
2460   //
2461   //  Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2462   //
2463
2464
2465   TObjArray *arr=AliTPCLaserTrack::GetTracks();
2466   TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2467   TVector3 vDir;
2468   TVector3 vSt;
2469
2470   Int_t firstbeam=0;
2471   Int_t lastbeam=336/2;
2472   if ( (sector/18)%2 ) {
2473     firstbeam=336/2;
2474     lastbeam=336;
2475   }
2476
2477   mindist=1000000;
2478   Int_t itrackMin=-1;
2479   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2480     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2481 //     if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2482     vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2483     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2484     if (TMath::Abs(deltaZ)>40) continue;
2485     vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2486     vSt.RotateZ(ltr->GetAlpha());
2487     vDir.RotateZ(ltr->GetAlpha());
2488
2489     Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2490
2491     if (dist<mindist){
2492       mindist=dist;
2493       itrackMin=itrack;
2494     }
2495
2496   }
2497   itrackMin2=-1;
2498   Float_t mindist2=10;
2499   for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2500     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);  //get the track
2501     if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2502
2503     Double_t deltaZ=ltr->GetZ()-peakpos[2];
2504     if (TMath::Abs(deltaZ)>40) continue;
2505
2506     Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2507     if (dist>1) continue;
2508
2509     if (dist<mindist2){
2510       mindist2=dist;
2511       itrackMin2=itrack;
2512     }
2513   }
2514   mindist=mindist2;
2515   return itrackMin2;
2516
2517 }
2518
2519 //_____________________________________________________________________
2520 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2521 {
2522   //
2523   // return true if pad is on the edge of a row
2524   //
2525   Int_t edge1   = 0;
2526   if ( pad == edge1 ) return kTRUE;
2527   Int_t edge2   = fROC->GetNPads(sector,row)-1;
2528   if ( pad == edge2 ) return kTRUE;
2529
2530   return kFALSE;
2531 }
2532
2533 //_____________________________________________________________________
2534 TObjArray* AliTPCCalibCE::SetupMeasured()
2535 {
2536   //
2537   // setup array of measured laser tracks and CE
2538   //
2539
2540   TObjArray *arrIdeal    = AliTPCLaserTrack::GetTracks();
2541   TObjArray *arrMeasured = new TObjArray(338);
2542   arrMeasured->SetOwner();
2543   for(Int_t itrack=0;itrack<336;itrack++){
2544     arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2545   }
2546
2547   //"tracks" for CE
2548   AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2549   ltrce->SetId(336);
2550   ltrce->SetSide(0);
2551   ltrce->fVecSec=new TVectorD(557568/2);
2552   ltrce->fVecP2=new TVectorD(557568/2);
2553   ltrce->fVecPhi=new TVectorD(557568/2);
2554   ltrce->fVecGX=new TVectorD(557568/2);
2555   ltrce->fVecGY=new TVectorD(557568/2);
2556   ltrce->fVecGZ=new TVectorD(557568/2);
2557   ltrce->fVecLX=new TVectorD(557568/2);
2558   ltrce->fVecLY=new TVectorD(557568/2);
2559   ltrce->fVecLZ=new TVectorD(557568/2);
2560
2561   arrMeasured->AddAt(ltrce,336); //CE A-Side
2562
2563   ltrce=new AliTPCLaserTrack;
2564   ltrce->SetId(337);
2565   ltrce->SetSide(1);
2566   ltrce->fVecSec=new TVectorD(557568/2);
2567   ltrce->fVecP2=new TVectorD(557568/2);
2568   ltrce->fVecPhi=new TVectorD(557568/2);
2569   ltrce->fVecGX=new TVectorD(557568/2);
2570   ltrce->fVecGY=new TVectorD(557568/2);
2571   ltrce->fVecGZ=new TVectorD(557568/2);
2572   ltrce->fVecLX=new TVectorD(557568/2);
2573   ltrce->fVecLY=new TVectorD(557568/2);
2574   ltrce->fVecLZ=new TVectorD(557568/2);
2575   arrMeasured->AddAt(ltrce,337); //CE C-Side
2576
2577   return arrMeasured;
2578 }
2579
2580 //_____________________________________________________________________
2581 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2582 {
2583   //
2584   // reset array of measured laser tracks and CE
2585   //
2586   for(Int_t itrack=0;itrack<338;itrack++){
2587     AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2588     ltr->fVecSec->Zero();
2589     ltr->fVecP2->Zero();
2590     ltr->fVecPhi->Zero();
2591     ltr->fVecGX->Zero();
2592     ltr->fVecGY->Zero();
2593     ltr->fVecGZ->Zero();
2594     ltr->fVecLX->Zero();
2595     ltr->fVecLY->Zero();
2596     ltr->fVecLZ->Zero();
2597   }
2598 }
2599
2600 //_____________________________________________________________________
2601 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2602 {
2603   //
2604   // Add ideal CE positions to the ideal track data
2605   //
2606
2607   arr->Expand(338);
2608   //"tracks" for CE
2609   AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2610   ltrceA->SetId(336);
2611   ltrceA->SetSide(0);
2612   ltrceA->fVecSec=new TVectorD(557568/2);
2613   ltrceA->fVecP2=new TVectorD(557568/2);
2614   ltrceA->fVecPhi=new TVectorD(557568/2);
2615   ltrceA->fVecGX=new TVectorD(557568/2);
2616   ltrceA->fVecGY=new TVectorD(557568/2);
2617   ltrceA->fVecGZ=new TVectorD(557568/2);
2618   ltrceA->fVecLX=new TVectorD(557568/2);
2619   ltrceA->fVecLY=new TVectorD(557568/2);
2620   ltrceA->fVecLZ=new TVectorD(557568/2);
2621   arr->AddAt(ltrceA,336); //CE A-Side
2622
2623   AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2624   ltrceC->SetId(337);
2625   ltrceC->SetSide(1);
2626   ltrceC->fVecSec=new TVectorD(557568/2);
2627   ltrceC->fVecP2=new TVectorD(557568/2);
2628   ltrceC->fVecPhi=new TVectorD(557568/2);
2629   ltrceC->fVecGX=new TVectorD(557568/2);
2630   ltrceC->fVecGY=new TVectorD(557568/2);
2631   ltrceC->fVecGZ=new TVectorD(557568/2);
2632   ltrceC->fVecLX=new TVectorD(557568/2);
2633   ltrceC->fVecLY=new TVectorD(557568/2);
2634   ltrceC->fVecLZ=new TVectorD(557568/2);
2635   arr->AddAt(ltrceC,337); //CE C-Side
2636
2637   //Calculate ideal positoins
2638   Float_t gxyz[3];
2639   Float_t lxyz[3];
2640   Int_t channelSideA=0;
2641   Int_t channelSideC=0;
2642   Int_t channelSide=0;
2643   AliTPCLaserTrack *ltrce=0x0;
2644
2645   for (Int_t isector=0; isector<72; ++isector){
2646     Int_t side=((isector/18)%2);
2647     for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2648       for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2649         fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2650         fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2651         if (side==0) {
2652           ltrce=ltrceA;
2653           channelSide=channelSideA;
2654         } else {
2655           ltrce=ltrceC;
2656           channelSide=channelSideC;
2657         }
2658
2659         ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2660         ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2661         ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2662         ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2663         ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2664 //         ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2665         ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2666         ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2667 //         ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2668
2669         if (side==0){
2670           ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2671           ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2672           ++channelSideA;
2673         }
2674         else {
2675           ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2676           ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2677           ++channelSideC;
2678         }
2679       }
2680     }
2681   }
2682
2683
2684 }
2685
2686 //_____________________________________________________________________
2687 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2688 {
2689   //
2690   // calculate the drift velocity from the reconstructed clusters associated
2691   // to the ideal laser tracks
2692   // use two different fit scenarios: Separate fit for A- and C-Side
2693   //                                  Common fit for A- and C-Side
2694   //
2695
2696   if (!fArrFitGraphs){
2697     fArrFitGraphs=new TObjArray;
2698   }
2699
2700 //   static TLinearFitter fdriftA(5,"hyp4");
2701 //   static TLinearFitter fdriftC(5,"hyp4");
2702 //   static TLinearFitter fdriftAC(6,"hyp5");
2703   Double_t timestamp=fTimeBursts[burst];
2704
2705   static TLinearFitter fdriftA(4,"hyp3");
2706   static TLinearFitter fdriftC(4,"hyp3");
2707   static TLinearFitter fdriftAC(5,"hyp4");
2708   TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2709
2710   Float_t chi2A = 10;
2711   Float_t chi2C = 10;
2712   Float_t chi2AC = 10;
2713   Int_t npointsA=0;
2714   Int_t npointsC=0;
2715   Int_t npointsAC=0;
2716
2717   Double_t minres[3]={20.,1,0.8};
2718   //----
2719   for(Int_t i=0;i<3;i++){
2720
2721     fdriftA.ClearPoints();
2722     fdriftC.ClearPoints();
2723     fdriftAC.ClearPoints();
2724
2725     chi2A = 10;
2726     chi2C = 10;
2727     chi2AC = 10;
2728     npointsA=0;
2729     npointsC=0;
2730     npointsAC=0;
2731
2732     for (Int_t itrack=0; itrack<338; ++itrack){
2733       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2734       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2735
2736       //-- Exclude the tracks which has the biggest inclanation angle
2737       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2738       Int_t clustercounter=0;
2739       Int_t indexMax=159;
2740
2741       //-- exclude the low intensity tracks
2742
2743       for (Int_t index=0; index<indexMax; ++index){
2744
2745         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2746         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2747         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2748
2749         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2750       }
2751       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2752       clustercounter=0;
2753
2754
2755       //-- drift length
2756       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2757
2758       if (itrack>335) indexMax=557568/2;
2759       for (Int_t index=0; index<indexMax; ++index){
2760         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2761         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2762         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2763         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2764
2765         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2766         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2767         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2768         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2769
2770       //cut if no track info available
2771         if (iltr->GetBundle()==0) continue;
2772         if (iR<133||mR<133) continue;
2773         if(TMath::Abs(mltr->fVecP2->GetMatrixArray()[index])>minres[i]) continue;
2774
2775         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2776         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2777
2778         //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2779         Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2780
2781         if (iltr->GetSide()==0){
2782           fdriftA.AddPoint(xxx,mdrift,1);
2783         }else{
2784           fdriftC.AddPoint(xxx,mdrift,1);
2785         }
2786 //         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2787         Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2788         fdriftAC.AddPoint(xxx2,mdrift,1);
2789
2790       }//end index loop
2791     }//end laser track loop
2792
2793   //perform fit
2794     fdriftA.Eval();
2795     fdriftC.Eval();
2796     fdriftAC.Eval();
2797
2798
2799
2800   //get fit values
2801     fdriftA.GetParameters(fitA);
2802     fdriftC.GetParameters(fitC);
2803     fdriftAC.GetParameters(fitAC);
2804
2805   //Parameters:  0 linear offset
2806   //             1 mean drift velocity correction factor
2807   //             2 relative global y gradient
2808   //             3 radial deformation
2809   //             4 IROC/OROC offset
2810
2811 //      FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2812
2813     for (Int_t itrack=0; itrack<338; ++itrack){
2814       AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2815       AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2816
2817       //-- Exclude the tracks which has the biggest inclanation angle
2818       if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2819       Int_t clustercounter=0;
2820       Int_t indexMax=159;
2821
2822       //-- exclude the low intensity tracks
2823
2824       for (Int_t index=0; index<indexMax; ++index){
2825         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2826         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2827         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2828         if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2829       }
2830       if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2831       clustercounter=0;
2832
2833       //-- drift length
2834       Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2835
2836       if (itrack>335) indexMax=557568/2;
2837       for (Int_t index=0; index<indexMax; ++index){
2838         Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2839         Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2840         Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2841         Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2842
2843         Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2844         Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2845         Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2846         Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2847
2848       //cut if no track info available
2849         if (iR<60||mR<60) continue;
2850
2851         Double_t ldrift  = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2852         Double_t mdrift  = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2853
2854         TVectorD *v=&fitA;
2855         if (iltr->GetSide()==1) v=&fitC;
2856 //         Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR)+(*v)[4]*( iltr->fVecSec->GetMatrixArray()[index]>35);
2857         Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2858
2859         mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2860
2861       }
2862     }
2863
2864     fitA.ResizeTo(7);
2865     fitC.ResizeTo(7);
2866     fitAC.ResizeTo(8);
2867
2868 //set statistics values
2869
2870     npointsA= fdriftA.GetNpoints();
2871     if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2872     fitA[5]=npointsA;
2873     fitA[6]=chi2A;
2874
2875     npointsC= fdriftC.GetNpoints();
2876     if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2877     fitC[5]=npointsC;
2878     fitC[6]=chi2C;
2879
2880     npointsAC= fdriftAC.GetNpoints();
2881     if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2882     fitAC[5]=npointsAC;
2883     fitAC[6]=chi2AC;
2884
2885     if (fStreamLevel>2){
2886     //laser track information
2887       (*GetDebugStreamer()) << "DriftV" <<
2888         "iter="   << i <<
2889         "run="    << fRunNumber <<
2890         "time="   << timestamp <<
2891         "fitA.="  << &fitA <<
2892         "fitC.="  << &fitC <<
2893         "fitAC.=" << &fitAC <<
2894         "\n";
2895
2896
2897     }
2898
2899   }
2900 //-----
2901
2902
2903   //Parameters:  0 linear offset (global)
2904   //             1 mean drift velocity correction factor
2905   //             2 relative global y gradient
2906   //             3 radial deformation
2907   //             4 IROC/OROC offset
2908   //             5 linear offset relative A-C
2909
2910   //get graphs
2911   TGraphErrors *grA[7];
2912   TGraphErrors *grC[7];
2913   TGraphErrors *grAC[8];
2914   TString names("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_");
2915   TString namesAC("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_;GRAPH_MEAN_DELAYC_LASER_ALL_");
2916
2917   TObjArray *arrNames=names.Tokenize(";");
2918   TObjArray *arrNamesAC=namesAC.Tokenize(";");
2919
2920   //A-Side graphs
2921   for (Int_t i=0; i<7; ++i){
2922     TString grName=arrNames->UncheckedAt(i)->GetName();
2923     grName+="A";
2924     grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2925     if (!grA[i]){
2926       grA[i]=new TGraphErrors;
2927       grA[i]->SetName(grName.Data());
2928       grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2929       fArrFitGraphs->Add(grA[i]);
2930     }
2931 //     Int_t ipoint=grA[i]->GetN();
2932     Int_t ipoint=burst;
2933     grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2934     grA[i]->SetPointError(ipoint,60,.0001);
2935     if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2936   }
2937
2938   //C-Side graphs
2939   for (Int_t i=0; i<7; ++i){
2940     TString grName=arrNames->UncheckedAt(i)->GetName();
2941     grName+="C";
2942     grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2943     if (!grC[i]){
2944       grC[i]=new TGraphErrors;
2945       grC[i]->SetName(grName.Data());
2946       grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2947       fArrFitGraphs->Add(grC[i]);
2948     }
2949 //     Int_t ipoint=grC[i]->GetN();
2950     Int_t ipoint=burst;
2951     grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2952     grC[i]->SetPointError(ipoint,60,.0001);
2953     if (i<4) grC[i]->SetPointError(ipoint,60,fdriftC.GetCovarianceMatrixElement(i,i));
2954   }
2955
2956   //AC-Side graphs
2957   for (Int_t i=0; i<8; ++i){
2958     TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2959     grName+="AC";
2960     grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2961     if (!grAC[i]){
2962       grAC[i]=new TGraphErrors;
2963       grAC[i]->SetName(grName.Data());
2964       grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2965       fArrFitGraphs->Add(grAC[i]);
2966     }
2967 //     Int_t ipoint=grAC[i]->GetN();
2968     Int_t ipoint=burst;
2969     grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2970     grAC[i]->SetPointError(ipoint,60,.0001);
2971     if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftAC.GetCovarianceMatrixElement(i,i));
2972   }
2973
2974   if (fDebugLevel>10){
2975     printf("A side fit parameters:\n");
2976     fitA.Print();
2977     printf("\nC side fit parameters:\n");
2978     fitC.Print();
2979     printf("\nAC side fit parameters:\n");
2980     fitAC.Print();
2981   }
2982   delete arrNames;
2983   delete arrNamesAC;
2984 }
2985
2986 //_____________________________________________________________________
2987 Double_t AliTPCCalibCE::SetBurstHnDrift()
2988 {
2989   //
2990   // Create a new THnSparse for the current burst
2991   // return the time of the current burst
2992   //
2993   Int_t i=0;
2994   for(i=0; i<fTimeBursts.GetNrows(); ++i){
2995     if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2996     if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2997       fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2998       return fTimeBursts(i);
2999     }
3000   }
3001
3002   CreateDVhist();
3003   fArrHnDrift.AddAt(fHnDrift,i);
3004   fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
3005   for (i=0;i<14;++i){
3006     fPeaks[i]=0;
3007     fPeakWidths[i]=0;
3008   }
3009   fEventInBunch=0;
3010   return fTimeStamp;
3011 }
3012
3013 //_____________________________________________________________________
3014 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
3015 {
3016   //
3017   //  Write class to file
3018   //  option can be specified in the dir option:
3019   //  options:
3020   //    name=<objname>: the name of the calibration object in file will be <objname>
3021   //    type=<type>:    the saving type:
3022   //                    0 - write the complte object
3023   //                    1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
3024   //                    2 - like 2, but in addition delete objects that will most probably not be used for calibration
3025   //                    3 - store only calibration output, don't store the reference histograms
3026   //                        and THnSparse (requires Analyse called before)
3027   //
3028   //  NOTE: to read the object back, the ReadFromFile function should be used
3029   //
3030
3031   TString sDir(dir);
3032   TString objName=GetName();
3033   Int_t type=0;
3034
3035   //get options
3036   TObjArray *arr=sDir.Tokenize(",");
3037   TIter next(arr);
3038   TObjString *s;
3039   while ( (s=(TObjString*)next()) ){
3040     TString optString=s->GetString();
3041     optString.Remove(TString::kBoth,' ');
3042     if (optString.BeginsWith("name=")){
3043       objName=optString.ReplaceAll("name=","");
3044     }
3045     if (optString.BeginsWith("type=")){
3046       optString.ReplaceAll("type=","");
3047       type=optString.Atoi();
3048     }
3049   }
3050
3051   if ( type==4 ){
3052     // only for the new algorithm
3053     AliTPCCalibCE ce;
3054     ce.fArrFitGraphs=fArrFitGraphs;
3055     ce.fNevents=fNevents;
3056     ce.fTimeBursts.ResizeTo(fTimeBursts.GetNrows());
3057     ce.fTimeBursts=fTimeBursts;
3058     ce.fProcessNew=kTRUE;
3059     TFile f(filename,"recreate");
3060     ce.Write(objName.Data());
3061     fArrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3062     f.Close();
3063     ce.fArrFitGraphs=0x0;
3064     return;
3065   }
3066
3067
3068   if (type==1||type==2) {
3069     //delete most probably not needed stuff
3070     //This requires Analyse to be called after reading the object from file
3071     fCalRocArrayT0.Delete();
3072     fCalRocArrayT0Err.Delete();
3073     fCalRocArrayQ.Delete();
3074     fCalRocArrayRMS.Delete();
3075     fCalRocArrayOutliers.Delete();
3076   }
3077   if (type==2||type==3){
3078     fParamArrayEventPol1.Delete();
3079     fParamArrayEventPol2.Delete();
3080   }
3081
3082   TObjArray histoQArray(72);
3083   TObjArray histoT0Array(72);
3084   TObjArray histoRMSArray(72);
3085   TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3086
3087   //save all large 2D histograms in separte pointers
3088   //to have a smaller memory print when saving the object
3089   if (type==1||type==2||type==3){
3090     for (Int_t i=0; i<72; ++i){
3091       histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3092       histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3093       histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3094     }
3095     fHistoQArray.SetOwner(kFALSE);
3096     fHistoT0Array.SetOwner(kFALSE);
3097     fHistoRMSArray.SetOwner(kFALSE);
3098     fHistoQArray.Clear();
3099     fHistoT0Array.Clear();
3100     fHistoRMSArray.Clear();
3101
3102     for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3103       arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3104     }
3105     fArrHnDrift.SetOwner(kFALSE);
3106     fArrHnDrift.Clear();
3107   }
3108
3109
3110   TDirectory *backup = gDirectory;
3111
3112   TFile f(filename,"recreate");
3113   Write(objName.Data());
3114   if (type==1||type==2) {
3115     histoQArray.Write("histoQArray",TObject::kSingleKey);
3116     histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3117     histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3118     arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3119   }
3120
3121   f.Save();
3122   f.Close();
3123
3124   //move histograms back to the object
3125   if (type==1||type==2){
3126     for (Int_t i=0; i<72; ++i){
3127       fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3128       fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3129       fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3130     }
3131     fHistoQArray.SetOwner(kTRUE);
3132     fHistoT0Array.SetOwner(kTRUE);
3133     fHistoRMSArray.SetOwner(kTRUE);
3134
3135     for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3136       fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3137     }
3138     fArrHnDrift.SetOwner(kTRUE);
3139   }
3140
3141   if ( backup ) backup->cd();
3142 }
3143 //_____________________________________________________________________
3144 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3145 {
3146   //
3147   // Read object from file
3148   // Handle properly if the histogram arrays were stored separately
3149   // call Analyse to make sure to have the calibration relevant information in the object
3150   //
3151
3152   TFile f(filename);
3153   if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3154   TList *l=f.GetListOfKeys();
3155   TIter next(l);
3156   TKey *key=0x0;
3157   TObject *o=0x0;
3158
3159   AliTPCCalibCE *ce=0x0;
3160   TObjArray *histoQArray=0x0;
3161   TObjArray *histoT0Array=0x0;
3162   TObjArray *histoRMSArray=0x0;
3163   TObjArray *arrHnDrift=0x0;
3164
3165   while ( (key=(TKey*)next()) ){
3166     o=key->ReadObj();
3167     if ( o->IsA()==AliTPCCalibCE::Class() ){
3168       ce=(AliTPCCalibCE*)o;
3169     } else if ( o->IsA()==TObjArray::Class() ){
3170       TString name=key->GetName();
3171       if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3172       if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3173       if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3174       if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3175     }
3176   }
3177
3178   if (ce){
3179   //move histograms back to the object
3180     TH1* hist=0x0;
3181     if (histoQArray){
3182       for (Int_t i=0; i<72; ++i){
3183         hist=(TH1*)histoQArray->UncheckedAt(i);
3184         if (hist) hist->SetDirectory(0x0);
3185         ce->fHistoQArray.AddAt(hist,i);
3186       }
3187       ce->fHistoQArray.SetOwner(kTRUE);
3188     }
3189
3190     if (histoT0Array) {
3191       for (Int_t i=0; i<72; ++i){
3192         hist=(TH1*)histoT0Array->UncheckedAt(i);
3193         if (hist) hist->SetDirectory(0x0);
3194         ce->fHistoT0Array.AddAt(hist,i);
3195       }
3196       ce->fHistoT0Array.SetOwner(kTRUE);
3197     }
3198
3199     if (histoRMSArray){
3200       for (Int_t i=0; i<72; ++i){
3201         hist=(TH1*)histoRMSArray->UncheckedAt(i);
3202         if (hist) hist->SetDirectory(0x0);
3203         ce->fHistoRMSArray.AddAt(hist,i);
3204       }
3205       ce->fHistoRMSArray.SetOwner(kTRUE);
3206     }
3207
3208     if (arrHnDrift){
3209       for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3210         THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3211         ce->fArrHnDrift.AddAt(hSparse,i);
3212       }
3213     }
3214
3215     ce->Analyse();
3216   }
3217   f.Close();
3218
3219   return ce;
3220 }