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