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