#include <TObjArray.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2S.h>
#include <TF1.h>
#include <TString.h>
#include <TVectorF.h>
#include <TVectorD.h>
#include <TVector3.h>
#include <TMatrixD.h>
#include <TMath.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TString.h>
#include <TMap.h>
#include <TDirectory.h>
#include <TSystem.h>
#include <TFile.h>
#include <TCollection.h>
#include <TTimeStamp.h>
#include <TList.h>
#include <TKey.h>
#include <TSpectrum.h>
#include "AliLog.h"
#include "AliRawReader.h"
#include "AliRawReaderRoot.h"
#include "AliRawReaderDate.h"
#include "AliRawEventHeaderBase.h"
#include "AliTPCCalROC.h"
#include "AliTPCCalPad.h"
#include "AliTPCROC.h"
#include "AliTPCParam.h"
#include "AliTPCCalibCE.h"
#include "AliMathBase.h"
#include "AliTPCTransform.h"
#include "AliTPCLaserTrack.h"
#include "TTreeStream.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "event.h"
ClassImp(AliTPCCalibCE)
AliTPCCalibCE::AliTPCCalibCE() :
AliTPCCalibRawBase(),
fNbinsT0(200),
fXminT0(-5),
fXmaxT0(5),
fNbinsQ(200),
fXminQ(1),
fXmaxQ(40),
fNbinsRMS(100),
fXminRMS(0.1),
fXmaxRMS(5.1),
fPeakDetMinus(2),
fPeakDetPlus(3),
fPeakIntMinus(2),
fPeakIntPlus(2),
fNoiseThresholdMax(5.),
fNoiseThresholdSum(8.),
fIsZeroSuppressed(kFALSE),
fLastSector(-1),
fSecRejectRatio(.4),
fParam(new AliTPCParam),
fPedestalTPC(0x0),
fPadNoiseTPC(0x0),
fPedestalROC(0x0),
fPadNoiseROC(0x0),
fCalRocArrayT0(72),
fCalRocArrayT0Err(72),
fCalRocArrayQ(72),
fCalRocArrayRMS(72),
fCalRocArrayOutliers(72),
fHistoQArray(72),
fHistoT0Array(72),
fHistoRMSArray(72),
fMeanT0rms(0),
fMeanQrms(0),
fMeanRMSrms(0),
fHistoTmean(72),
fParamArrayEventPol1(72),
fParamArrayEventPol2(72),
fTMeanArrayEvent(72),
fQMeanArrayEvent(72),
fVEventTime(1000),
fVEventNumber(1000),
fVTime0SideA(1000),
fVTime0SideC(1000),
fEventId(-1),
fOldRunNumber(0),
fPadTimesArrayEvent(72),
fPadQArrayEvent(72),
fPadRMSArrayEvent(72),
fPadPedestalArrayEvent(72),
fCurrentChannel(-1),
fCurrentSector(-1),
fCurrentRow(-1),
fMaxPadSignal(-1),
fMaxTimeBin(-1),
fPadPedestal(0),
fPadNoise(0),
fVTime0Offset(72),
fVTime0OffsetCounter(72),
fVMeanQ(72),
fVMeanQCounter(72),
fCurrentCETimeRef(0),
fProcessOld(kTRUE),
fProcessNew(kFALSE),
fAnalyseNew(kTRUE),
fHnDrift(0x0),
fArrHnDrift(100),
fTimeBursts(100),
fArrFitGraphs(0x0),
fEventInBunch(0)
{
SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
fFirstTimeBin=650;
fLastTimeBin=1030;
fParam->Update();
for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
for (Int_t i=0;i<14;++i){
fPeaks[i]=0;
fPeakWidths[i]=0;
}
for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
}
AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
AliTPCCalibRawBase(sig),
fNbinsT0(sig.fNbinsT0),
fXminT0(sig.fXminT0),
fXmaxT0(sig.fXmaxT0),
fNbinsQ(sig.fNbinsQ),
fXminQ(sig.fXminQ),
fXmaxQ(sig.fXmaxQ),
fNbinsRMS(sig.fNbinsRMS),
fXminRMS(sig.fXminRMS),
fXmaxRMS(sig.fXmaxRMS),
fPeakDetMinus(sig.fPeakDetMinus),
fPeakDetPlus(sig.fPeakDetPlus),
fPeakIntMinus(sig.fPeakIntMinus),
fPeakIntPlus(sig.fPeakIntPlus),
fNoiseThresholdMax(sig.fNoiseThresholdMax),
fNoiseThresholdSum(sig.fNoiseThresholdSum),
fIsZeroSuppressed(sig.fIsZeroSuppressed),
fLastSector(-1),
fSecRejectRatio(.4),
fParam(new AliTPCParam),
fPedestalTPC(0x0),
fPadNoiseTPC(0x0),
fPedestalROC(0x0),
fPadNoiseROC(0x0),
fCalRocArrayT0(72),
fCalRocArrayT0Err(72),
fCalRocArrayQ(72),
fCalRocArrayRMS(72),
fCalRocArrayOutliers(72),
fHistoQArray(72),
fHistoT0Array(72),
fHistoRMSArray(72),
fMeanT0rms(sig.fMeanT0rms),
fMeanQrms(sig.fMeanQrms),
fMeanRMSrms(sig.fMeanRMSrms),
fHistoTmean(72),
fParamArrayEventPol1(72),
fParamArrayEventPol2(72),
fTMeanArrayEvent(72),
fQMeanArrayEvent(72),
fVEventTime(sig.fVEventTime),
fVEventNumber(sig.fVEventNumber),
fVTime0SideA(sig.fVTime0SideA),
fVTime0SideC(sig.fVTime0SideC),
fEventId(-1),
fOldRunNumber(0),
fPadTimesArrayEvent(72),
fPadQArrayEvent(72),
fPadRMSArrayEvent(72),
fPadPedestalArrayEvent(72),
fCurrentChannel(-1),
fCurrentSector(-1),
fCurrentRow(-1),
fMaxPadSignal(-1),
fMaxTimeBin(-1),
fPadPedestal(0),
fPadNoise(0),
fVTime0Offset(72),
fVTime0OffsetCounter(72),
fVMeanQ(72),
fVMeanQCounter(72),
fCurrentCETimeRef(0),
fProcessOld(sig.fProcessOld),
fProcessNew(sig.fProcessNew),
fAnalyseNew(sig.fAnalyseNew),
fHnDrift(0x0),
fArrHnDrift(100),
fTimeBursts(100),
fArrFitGraphs(0x0),
fEventInBunch(0)
{
for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
for (Int_t iSec = 0; iSec < 72; ++iSec){
const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
if ( hQ != 0x0 ){
TH2S *hNew = new TH2S(*hQ);
hNew->SetDirectory(0);
fHistoQArray.AddAt(hNew,iSec);
}
if ( hT0 != 0x0 ){
TH2S *hNew = new TH2S(*hT0);
hNew->SetDirectory(0);
fHistoT0Array.AddAt(hNew,iSec);
}
if ( hRMS != 0x0 ){
TH2S *hNew = new TH2S(*hRMS);
hNew->SetDirectory(0);
fHistoRMSArray.AddAt(hNew,iSec);
}
}
TObjArray *arr=0x0;
for (Int_t iSec=0; iSec<72; ++iSec){
arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
if ( arr ){
TObjArray *arrEvents = new TObjArray(arr->GetSize());
fParamArrayEventPol1.AddAt(arrEvents, iSec);
for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
arrEvents->AddAt(new TVectorD(*vec),iEvent);
}
arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
if ( arr ){
TObjArray *arrEvents = new TObjArray(arr->GetSize());
fParamArrayEventPol2.AddAt(arrEvents, iSec);
for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
arrEvents->AddAt(new TVectorD(*vec),iEvent);
}
TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
if ( vMeanTime )
fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
if ( vMeanQ )
fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
}
fVEventTime.ResizeTo(sig.fVEventTime);
fVEventNumber.ResizeTo(sig.fVEventNumber);
fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
fParam->Update();
for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
TObject *o=sig.fArrHnDrift.UncheckedAt(i);
if (o){
TObject *newo=o->Clone("fHnDrift");
fArrHnDrift.AddAt(newo,i);
if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
}
}
for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
fTimeBursts[i]=sig.fTimeBursts[i];
}
for (Int_t i=0;i<14;++i){
fPeaks[i]=sig.fPeaks[i];
fPeakWidths[i]=sig.fPeakWidths[i];
}
if (sig.fArrFitGraphs) {
fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
fArrFitGraphs->SetOwner();
}
for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
}
AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
AliTPCCalibRawBase(),
fNbinsT0(200),
fXminT0(-5),
fXmaxT0(5),
fNbinsQ(200),
fXminQ(1),
fXmaxQ(40),
fNbinsRMS(100),
fXminRMS(0.1),
fXmaxRMS(5.1),
fPeakDetMinus(2),
fPeakDetPlus(3),
fPeakIntMinus(2),
fPeakIntPlus(2),
fNoiseThresholdMax(5.),
fNoiseThresholdSum(8.),
fIsZeroSuppressed(kFALSE),
fLastSector(-1),
fSecRejectRatio(.4),
fParam(new AliTPCParam),
fPedestalTPC(0x0),
fPadNoiseTPC(0x0),
fPedestalROC(0x0),
fPadNoiseROC(0x0),
fCalRocArrayT0(72),
fCalRocArrayT0Err(72),
fCalRocArrayQ(72),
fCalRocArrayRMS(72),
fCalRocArrayOutliers(72),
fHistoQArray(72),
fHistoT0Array(72),
fHistoRMSArray(72),
fMeanT0rms(0),
fMeanQrms(0),
fMeanRMSrms(0),
fHistoTmean(72),
fParamArrayEventPol1(72),
fParamArrayEventPol2(72),
fTMeanArrayEvent(72),
fQMeanArrayEvent(72),
fVEventTime(1000),
fVEventNumber(1000),
fVTime0SideA(1000),
fVTime0SideC(1000),
fEventId(-1),
fOldRunNumber(0),
fPadTimesArrayEvent(72),
fPadQArrayEvent(72),
fPadRMSArrayEvent(72),
fPadPedestalArrayEvent(72),
fCurrentChannel(-1),
fCurrentSector(-1),
fCurrentRow(-1),
fMaxPadSignal(-1),
fMaxTimeBin(-1),
fPadPedestal(0),
fPadNoise(0),
fVTime0Offset(72),
fVTime0OffsetCounter(72),
fVMeanQ(72),
fVMeanQCounter(72),
fCurrentCETimeRef(0),
fProcessOld(kTRUE),
fProcessNew(kFALSE),
fAnalyseNew(kTRUE),
fHnDrift(0x0),
fArrHnDrift(100),
fTimeBursts(100),
fArrFitGraphs(0x0),
fEventInBunch(0)
{
SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
fFirstTimeBin=650;
fLastTimeBin=1030;
if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
for (Int_t i=0;i<14;++i){
fPeaks[i]=0;
fPeakWidths[i]=0;
}
fParam->Update();
for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
}
AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
{
if (&source == this) return *this;
new (this) AliTPCCalibCE(source);
return *this;
}
AliTPCCalibCE::~AliTPCCalibCE()
{
fCalRocArrayT0.Delete();
fCalRocArrayT0Err.Delete();
fCalRocArrayQ.Delete();
fCalRocArrayRMS.Delete();
fCalRocArrayOutliers.Delete();
fHistoQArray.Delete();
fHistoT0Array.Delete();
fHistoRMSArray.Delete();
fHistoTmean.Delete();
fParamArrayEventPol1.Delete();
fParamArrayEventPol2.Delete();
fTMeanArrayEvent.Delete();
fQMeanArrayEvent.Delete();
fPadTimesArrayEvent.Delete();
fPadQArrayEvent.Delete();
fPadRMSArrayEvent.Delete();
fPadPedestalArrayEvent.Delete();
fArrHnDrift.SetOwner();
fArrHnDrift.Delete();
if (fArrFitGraphs){
fArrFitGraphs->SetOwner();
delete fArrFitGraphs;
}
}
Int_t AliTPCCalibCE::Update(const Int_t icsector,
const Int_t icRow,
const Int_t icPad,
const Int_t icTimeBin,
const Float_t csignal)
{
if (!fProcessOld) return 0;
if (icRow<0) return 0;
if (icPad<0) return 0;
if (icTimeBin<0) return 0;
if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad;
if ( fCurrentChannel == -1 ) {
fLastSector=-1;
fCurrentChannel = iChannel;
fCurrentSector = icsector;
fCurrentRow = icRow;
}
if ( iChannel != fCurrentChannel ){
ProcessPad();
fLastSector=fCurrentSector;
fCurrentChannel = iChannel;
fCurrentSector = icsector;
fCurrentRow = icRow;
}
fPadSignal[icTimeBin]=csignal;
if ( csignal > fMaxPadSignal ){
fMaxPadSignal = csignal;
fMaxTimeBin = icTimeBin;
}
return 0;
}
void AliTPCCalibCE::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)
{
if (!fProcessNew) return;
if (sector<36) return;
if (row<40) return;
if (length<3||length>10) return;
UShort_t timeBin = (UShort_t)startTimeBin;
if (timeBin<280) return;
Double_t timeBurst=SetBurstHnDrift();
Int_t cePeak=((sector/18)%2)*7+6;
if (fEventInBunch==1 && fPeaks[cePeak]==0) {
fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
FindLaserLayers();
fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
fHnDrift->Reset();
}
Int_t padFill=pad;
if (fEventInBunch==0||(fPeaks[cePeak]>100&&TMath::Abs((Short_t)fPeaks[cePeak]-(Short_t)timeBin)<(Short_t)fPeakWidths[cePeak])){
Int_t mod=5;
Int_t n=pad/mod;
padFill=mod*n+mod/2;
}
if (!IsPeakInRange(timeBin+length/2,sector)) return;
Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
(Double_t)padFill,(Double_t)timeBin,timeBurst};
for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
Float_t sig=(Float_t)signal[iTimeBin];
x[3]=timeBin;
fHnDrift->Fill(x,sig);
--timeBin;
}
}
void AliTPCCalibCE::FindLaserLayers()
{
for (Int_t iside=0;iside<2;++iside){
Int_t add=7*iside;
fHnDrift->GetAxis(0)->SetRangeUser(36+iside*18,53+iside*18);
TH1D *hproj=fHnDrift->Projection(3);
hproj->GetXaxis()->SetRangeUser(700,1030);
Int_t maxbin=hproj->GetMaximumBin();
Double_t binc=hproj->GetBinCenter(maxbin);
hproj->GetXaxis()->SetRangeUser(binc-5,binc+5);
fPeaks[add+6]=(UShort_t)TMath::Nint(binc);
fPeakWidths[add+6]=7;
hproj->GetXaxis()->SetRangeUser(0,maxbin-10);
TSpectrum s(6);
s.Search(hproj,2,"goff");
Int_t index[6];
TMath::Sort(6,s.GetPositionX(),index,kFALSE);
for (Int_t i=0; i<6; ++i){
fPeaks[i+add]=(UShort_t)TMath::Nint(s.GetPositionX()[index[i]]);
fPeakWidths[i+add]=5;
}
delete hproj;
}
}
void AliTPCCalibCE::FindPedestal(Float_t part)
{
Bool_t noPedestal = kTRUE;
if (fPedestalTPC&&fPadNoiseTPC){
if ( fCurrentSector!=fLastSector ){
fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
}
if ( fPedestalROC&&fPadNoiseROC ){
fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
noPedestal = kFALSE;
}
}
if ( noPedestal ) {
fPadPedestal = 0;
fPadNoise = 0;
if ( fIsZeroSuppressed ) return;
const Int_t kPedMax = 100;
Float_t max = 0;
Float_t maxPos = 0;
Int_t median = -1;
Int_t count0 = 0;
Int_t count1 = 0;
Float_t padSignal=0;
UShort_t histo[kPedMax];
memset(histo,0,kPedMax*sizeof(UShort_t));
for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
padSignal = fPadSignal[i];
if (padSignal<=0) continue;
if (padSignal>max && i>10) {
max = padSignal;
maxPos = i;
}
if (padSignal>kPedMax-1) continue;
histo[int(padSignal+0.5)]++;
count0++;
}
for (Int_t i=1; i<kPedMax; ++i){
if (count1<count0*0.5) median=i;
count1+=histo[i];
}
Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
for (Int_t idelta=1; idelta<10; ++idelta){
if (median-idelta<=0) continue;
if (median+idelta>kPedMax) continue;
if (count<part*count1){
count+=histo[median-idelta];
mean +=histo[median-idelta]*(median-idelta);
rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
count+=histo[median+idelta];
mean +=histo[median+idelta]*(median+idelta);
rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
}
}
if ( count > 0 ) {
mean/=count;
rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
fPadPedestal = mean;
fPadNoise = rms;
}
}
}
void AliTPCCalibCE::UpdateCETimeRef()
{
if ( fLastSector == fCurrentSector ) return;
Int_t sector=fCurrentSector;
if ( sector < 18 ) sector+=36;
fCurrentCETimeRef=0;
TVectorF *vtRef = GetTMeanEvents(sector);
if ( !vtRef ) return;
Int_t vtRefSize= vtRef->GetNrows();
if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
else vtRefSize=fNevents;
while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
fCurrentCETimeRef=(*vtRef)[vtRefSize];
AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
}
void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
{
Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
Int_t cemaxpos = 0;
Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise;
const Int_t kCemin = fPeakIntMinus;
const Int_t kCemax = fPeakIntPlus;
Float_t minDist = 25;
for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
Float_t tmean = fCurrentCETimeRef;
if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
minDist = tmean-maxima[imax];
cemaxpos = (Int_t)maxima[imax];
}
}
if (cemaxpos!=0){
ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
Float_t signal = fPadSignal[i]-fPadPedestal;
if (signal>0) {
ceTime+=signal*(i+0.5);
ceRMS +=signal*(i+0.5)*(i+0.5);
ceQsum+=signal;
}
}
}
}
if (ceQmax&&ceQsum>ceSumThreshold) {
ceTime/=ceQsum;
ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
ceTime-=GetL1PhaseTB();
fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;
fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
ceQsum/=norm;
fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
} else {
ceQmax=0;
ceTime=0;
ceRMS =0;
ceQsum=0;
}
param[0] = ceQmax;
param[1] = ceTime;
param[2] = ceRMS;
qSum = ceQsum;
}
Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
{
if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) )
++iTime2;
if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
}
return kTRUE;
}
void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
{
Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.));
Int_t count = 0;
for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
if (count<maxima.GetNrows()){
maxima.GetMatrixArray()[count++]=i;
GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
i-=(fPeakDetMinus+fPeakDetPlus-1);
}
}
}
}
void AliTPCCalibCE::ProcessPad()
{
FindPedestal();
TVectorF maxima(15);
FindLocalMaxima(maxima);
if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;
UpdateCETimeRef();
if ( fCurrentCETimeRef<1e-30 ) return;
TVectorD param(3);
Float_t qSum;
FindCESignal(param, qSum, maxima);
Double_t meanT = param[1];
Double_t sigmaT = param[2];
(*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
if ( GetStreamLevel()>0 ){
(*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
(*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
(*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
}
ResetPad();
}
void AliTPCCalibCE::EndEvent()
{
if (!fProcessOld) {
if (fProcessNew){
++fNevents;
++fEventInBunch;
}
return;
}
if ( fMaxTimeBin>-1 ) ProcessPad();
AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
TVectorD param(3);
TMatrixD dummy(3,3);
AliTPCCalROC *calIroc=new AliTPCCalROC(0);
AliTPCCalROC *calOroc=new AliTPCCalROC(36);
Double_t time0Side[2];
Double_t time0SideCount[2];
time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
for ( Int_t iSec = 36; iSec<72; ++iSec ){
time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
}
if ( time0SideCount[0] >0 )
time0Side[0]/=time0SideCount[0];
if ( time0SideCount[1] >0 )
time0Side[1]/=time0SideCount[1];
AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
Int_t nSecMeanT=0;
for ( Int_t iSec = 0; iSec<72; ++iSec ){
AliDebug(4,Form("Processing sector '%02d'\n",iSec));
TH1S *hMeanT = GetHistoTmean(iSec);
if ( !hMeanT ) continue;
if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
hMeanT->Reset();
AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
continue;
}
Double_t entries = hMeanT->GetEffectiveEntries();
Double_t sum = 0;
Short_t *arr = hMeanT->GetArray()+1;
Int_t ibin=0;
for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
sum+=arr[ibin];
if ( sum>=(entries/2.) ) break;
}
Int_t delta = 4;
Int_t firstBin = fFirstTimeBin+ibin-delta;
Int_t lastBin = fFirstTimeBin+ibin+delta;
if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
Int_t vSize=vMeanTime->GetNrows();
if ( vSize < fNevents+1 ){
vMeanTime->ResizeTo(vSize+100);
}
vSize=fVTime0SideA.GetNrows();
if ( vSize < fNevents+1 ){
fVTime0SideA.ResizeTo(vSize+100);
fVTime0SideC.ResizeTo(vSize+100);
}
fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
vMeanTime->GetMatrixArray()[fNevents]=median;
nSecMeanT++;
TVectorF *vTimes = GetPadTimesEvent(iSec);
if ( !vTimes ) continue;
AliTPCCalROC calIrocOutliers(0);
AliTPCCalROC calOrocOutliers(36);
TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
vSize=vMeanQ->GetNrows();
if ( vSize < fNevents+1 ){
vMeanQ->ResizeTo(vSize+100);
}
Float_t meanQ = 0;
if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
vMeanQ->GetMatrixArray()[fNevents]=meanQ;
for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
Float_t time = (*vTimes).GetMatrixArray()[iChannel];
if ( iSec < 36 ) {
calIroc->SetValue(iChannel, time);
if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
} else {
calOroc->SetValue(iChannel, time);
if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
}
if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
if ( GetStreamLevel()>0 ){
TTreeSRedirector *streamer=GetDebugStreamer();
if (streamer){
Int_t row=0;
Int_t pad=0;
Int_t padc=0;
Float_t q = (*GetPadQEvent(iSec))[iChannel];
Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
UInt_t channel=iChannel;
Int_t sector=iSec;
while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
pad = channel-fROC->GetRowIndexes(sector)[row];
padc = pad-(fROC->GetNPads(sector,row)/2);
Double_t t0Sec = 0;
if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
Double_t t0Side = time0Side[(iSec/18)%2];
(*streamer) << "DataPad" <<
"Event=" << fNevents <<
"RunNumber=" << fRunNumber <<
"TimeStamp=" << fTimeStamp <<
"Sector="<< sector <<
"Row=" << row<<
"Pad=" << pad <<
"PadC=" << padc <<
"PadSec="<< channel <<
"Time0Sec=" << t0Sec <<
"Time0Side=" << t0Side <<
"Time=" << time <<
"RMS=" << rms <<
"Sum=" << q <<
"MeanQ=" << meanQ <<
"\n";
}
}
}
if (GetDebugLevel()>0){
TVectorD paramPol1(3);
TVectorD paramPol2(6);
TMatrixD matPol1(3,3);
TMatrixD matPol2(6,6);
Float_t chi2Pol1=0;
Float_t chi2Pol2=0;
if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
if ( iSec < 36 ){
calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
} else {
calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
}
GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
}
if ( GetStreamLevel()>0 ){
TTreeSRedirector *streamer=GetDebugStreamer();
if ( streamer ) {
(*streamer) << "DataRoc" <<
"RunNumber=" << fRunNumber <<
"TimeStamp=" << fTimeStamp <<
"Sector="<< iSec <<
"hMeanT.=" << hMeanT <<
"median=" << median <<
"paramPol1.=" << ¶mPol1 <<
"paramPol2.=" << ¶mPol2 <<
"matPol1.=" << &matPol1 <<
"matPol2.=" << &matPol2 <<
"chi2Pol1=" << chi2Pol1 <<
"chi2Pol2=" << chi2Pol2 <<
"\n";
}
}
}
hMeanT->Reset();
}
if ( nSecMeanT == 0 ) return;
if ( fVEventTime.GetNrows() < fNevents+1 ) {
fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
}
fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
++fNevents;
if (fProcessNew) ++fEventInBunch;
fOldRunNumber = fRunNumber;
delete calIroc;
delete calOroc;
AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
}
TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
Int_t nbinsY, Float_t ymin, Float_t ymax,
const Char_t *type, Bool_t force)
{
if ( !force || arr->UncheckedAt(sector) )
return (TH2S*)arr->UncheckedAt(sector);
TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
nbinsY, ymin, ymax,
fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
hist->SetDirectory(0);
arr->AddAt(hist,sector);
return hist;
}
TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
{
TObjArray *arr = &fHistoT0Array;
return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
}
TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
{
TObjArray *arr = &fHistoQArray;
return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
}
TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
{
TObjArray *arr = &fHistoRMSArray;
return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
}
TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
const Char_t *type, Bool_t force)
{
if ( !force || arr->UncheckedAt(sector) )
return (TH1S*)arr->UncheckedAt(sector);
TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
hist->SetDirectory(0);
arr->AddAt(hist,sector);
return hist;
}
TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
{
TObjArray *arr = &fHistoTmean;
return GetHisto(sector, arr, "LastTmean", force);
}
TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
{
if ( !force || arr->UncheckedAt(sector) )
return (TVectorF*)arr->UncheckedAt(sector);
TVectorF *vect = new TVectorF(size);
arr->AddAt(vect,sector);
return vect;
}
TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
{
TObjArray *arr = &fPadTimesArrayEvent;
return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
}
TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
{
TObjArray *arr = &fPadQArrayEvent;
return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
}
TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
{
TObjArray *arr = &fPadRMSArrayEvent;
return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
}
TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
{
TObjArray *arr = &fPadPedestalArrayEvent;
return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
}
TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
{
TObjArray *arr = &fTMeanArrayEvent;
return GetVectSector(sector,arr,100,force);
}
TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
{
TObjArray *arr = &fQMeanArrayEvent;
return GetVectSector(sector,arr,100,force);
}
AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
{
if ( !force || arr->UncheckedAt(sector) )
return (AliTPCCalROC*)arr->UncheckedAt(sector);
AliTPCCalROC *croc = new AliTPCCalROC(sector);
arr->AddAt(croc,sector);
return croc;
}
AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
{
TObjArray *arr = &fCalRocArrayT0;
return GetCalRoc(sector, arr, force);
}
AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
{
TObjArray *arr = &fCalRocArrayT0Err;
return GetCalRoc(sector, arr, force);
}
AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
{
TObjArray *arr = &fCalRocArrayQ;
return GetCalRoc(sector, arr, force);
}
AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
{
TObjArray *arr = &fCalRocArrayRMS;
return GetCalRoc(sector, arr, force);
}
AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
{
TObjArray *arr = &fCalRocArrayOutliers;
return GetCalRoc(sector, arr, force);
}
TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
{
if ( !force || arr->UncheckedAt(sector) )
return (TObjArray*)arr->UncheckedAt(sector);
TObjArray *newArr = new TObjArray;
arr->AddAt(newArr,sector);
return newArr;
}
TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
{
TObjArray *arr = &fParamArrayEventPol1;
return GetParamArray(sector, arr, force);
}
TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
{
TObjArray *arr = &fParamArrayEventPol2;
return GetParamArray(sector, arr, force);
}
void AliTPCCalibCE::CreateDVhist()
{
TTimeStamp begin(2010,01,01,0,0,0);
TTimeStamp end(2035,01,01,0,0,0);
Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60;
Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
fHnDrift->Reset();
}
void AliTPCCalibCE::ResetEvent()
{
fLastSector=-1;
fCurrentSector=-1;
fCurrentRow=-1;
fCurrentChannel=-1;
ResetPad();
fPadTimesArrayEvent.Delete();
fPadQArrayEvent.Delete();
fPadRMSArrayEvent.Delete();
fPadPedestalArrayEvent.Delete();
for ( Int_t i=0; i<72; ++i ){
fVTime0Offset.GetMatrixArray()[i]=0;
fVTime0OffsetCounter.GetMatrixArray()[i]=0;
fVMeanQ.GetMatrixArray()[i]=0;
fVMeanQCounter.GetMatrixArray()[i]=0;
}
}
void AliTPCCalibCE::ResetPad()
{
for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
fPadSignal[i] = 0;
fMaxTimeBin = -1;
fMaxPadSignal = -1;
fPadPedestal = -1;
fPadNoise = -1;
}
void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
{
MergeBase(ce);
Int_t nCEevents = ce->GetNeventsProcessed();
if (fProcessOld&&ce->fProcessOld){
for (Int_t iSec=0; iSec<72; ++iSec){
TH2S *hRefQmerge = ce->GetHistoQ(iSec);
TH2S *hRefT0merge = ce->GetHistoT0(iSec);
TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
if ( hRefQmerge ){
TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
TH2S *hRefQ = GetHistoQ(iSec);
if ( hRefQ ) hRefQ->Add(hRefQmerge);
else {
TH2S *hist = new TH2S(*hRefQmerge);
hist->SetDirectory(0);
fHistoQArray.AddAt(hist, iSec);
}
hRefQmerge->SetDirectory(dir);
}
if ( hRefT0merge ){
TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
TH2S *hRefT0 = GetHistoT0(iSec);
if ( hRefT0 ) hRefT0->Add(hRefT0merge);
else {
TH2S *hist = new TH2S(*hRefT0merge);
hist->SetDirectory(0);
fHistoT0Array.AddAt(hist, iSec);
}
hRefT0merge->SetDirectory(dir);
}
if ( hRefRMSmerge ){
TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
TH2S *hRefRMS = GetHistoRMS(iSec);
if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
else {
TH2S *hist = new TH2S(*hRefRMSmerge);
hist->SetDirectory(0);
fHistoRMSArray.AddAt(hist, iSec);
}
hRefRMSmerge->SetDirectory(dir);
}
}
for (Int_t iSec=0; iSec<72; ++iSec){
TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
TObjArray *arrPol1 = 0x0;
TObjArray *arrPol2 = 0x0;
TVectorF *vMeanTime = 0x0;
TVectorF *vMeanQ = 0x0;
if ( arrPol1CE && arrPol2CE ){
arrPol1 = GetParamArrayPol1(iSec,kTRUE);
arrPol2 = GetParamArrayPol2(iSec,kTRUE);
arrPol1->Expand(fNevents+nCEevents);
arrPol2->Expand(fNevents+nCEevents);
}
if ( vMeanTimeCE && vMeanQCE ){
vMeanTime = GetTMeanEvents(iSec,kTRUE);
vMeanQ = GetQMeanEvents(iSec,kTRUE);
vMeanTime->ResizeTo(fNevents+nCEevents);
vMeanQ->ResizeTo(fNevents+nCEevents);
}
for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
if ( arrPol1CE && arrPol2CE ){
TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
if ( paramPol1 && paramPol2 ){
GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
}
}
if ( vMeanTimeCE && vMeanQCE ){
vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
}
}
}
const TVectorD& eventTimes = ce->fVEventTime;
const TVectorD& eventIds = ce->fVEventNumber;
const TVectorF& time0SideA = ce->fVTime0SideA;
const TVectorF& time0SideC = ce->fVTime0SideC;
fVEventTime.ResizeTo(fNevents+nCEevents);
fVEventNumber.ResizeTo(fNevents+nCEevents);
fVTime0SideA.ResizeTo(fNevents+nCEevents);
fVTime0SideC.ResizeTo(fNevents+nCEevents);
for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
Double_t evId = eventIds.GetMatrixArray()[iEvent];
Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
}
}
if (fProcessNew&&ce->fProcessNew) {
if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
AliError("Number of bursts in the instances to merge are different. No merging done!");
} else {
for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
if (h && hce) h->Add(hce);
else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
}
}
}
fNevents+=nCEevents;
}
Long64_t AliTPCCalibCE::Merge(TCollection * const list)
{
Long64_t nmerged=1;
TIter next(list);
AliTPCCalibCE *ce=0;
TObject *o=0;
while ( (o=next()) ){
ce=dynamic_cast<AliTPCCalibCE*>(o);
if (ce){
Merge(ce);
++nmerged;
}
}
return nmerged;
}
TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
{
TVectorD *xVar = 0x0;
TObjArray *aType = 0x0;
Int_t npoints=0;
if ( (sector<-2) || (sector>71) ) return 0x0;
if ( (xVariable<0) || (xVariable>2) ) return 0x0;
if ( (fitType<0) || (fitType>3) ) return 0x0;
if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0;
if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0;
if ( sector<0 && fitType!=2) return 0x0;
if (sector>=0){
if ( fitType==0 ){
if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
aType = &fParamArrayEventPol1;
if ( aType->At(sector)==0x0 ) return 0x0;
}
else if ( fitType==1 ){
if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
aType = &fParamArrayEventPol2;
if ( aType->At(sector)==0x0 ) return 0x0;
}
}
if ( xVariable == 0 ) xVar = &fVEventTime;
if ( xVariable == 1 ) xVar = &fVEventNumber;
if ( xVariable == 2 ) {
xVar = new TVectorD(fNevents);
for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
}
Double_t *x = new Double_t[fNevents];
Double_t *y = new Double_t[fNevents];
for (Int_t ievent =0; ievent<fNevents; ++ievent){
if ( fitType<2 ){
TObjArray *events = (TObjArray*)(aType->At(sector));
if ( events->GetSize()<=ievent ) break;
TVectorD *v = (TVectorD*)(events->At(ievent));
if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
} else if (fitType == 2) {
Double_t xValue=(*xVar)[ievent];
Double_t yValue=0;
if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
else if (sector==-1) yValue=fVTime0SideA(ievent);
else if (sector==-2) yValue=fVTime0SideC(ievent);
if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
}else if (fitType == 3) {
Double_t xValue=(*xVar)[ievent];
Double_t yValue=(*GetQMeanEvents(sector))[ievent];
if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
}
}
TGraph *gr = new TGraph(npoints);
Int_t *sortIndex = new Int_t[npoints];
TMath::Sort(npoints,x,sortIndex, kFALSE);
for (Int_t i=0;i<npoints;++i){
gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
}
if ( xVariable == 2 ) delete xVar;
delete [] x;
delete [] y;
delete [] sortIndex;
return gr;
}
void AliTPCCalibCE::Analyse()
{
if (fProcessOld){
TVectorD paramQ(3);
TVectorD paramT0(3);
TVectorD paramRMS(3);
TMatrixD dummy(3,3);
Float_t channelCounter=0;
fMeanT0rms=0;
fMeanQrms=0;
fMeanRMSrms=0;
for (Int_t iSec=0; iSec<72; ++iSec){
TH2S *hT0 = GetHistoT0(iSec);
if (!hT0 ) continue;
AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
TH2S *hQ = GetHistoQ(iSec);
TH2S *hRMS = GetHistoRMS(iSec);
Short_t *arrayhQ = hQ->GetArray();
Short_t *arrayhT0 = hT0->GetArray();
Short_t *arrayhRMS = hRMS->GetArray();
UInt_t nChannels = fROC->GetNChannels(iSec);
Int_t row=0;
Int_t pad=0;
Int_t padc=0;
for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
Float_t cogTime0 = -1000;
Float_t cogQ = -1000;
Float_t cogRMS = -1000;
Float_t cogOut = 0;
Float_t rms = 0;
Float_t rmsT0 = 0;
Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
fMeanQrms+=rms;
cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
fMeanT0rms+=rmsT0;
cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
fMeanRMSrms+=rms;
channelCounter++;
rocQ->SetValue(iChannel, cogQ*cogQ);
rocT0->SetValue(iChannel, cogTime0);
rocT0Err->SetValue(iChannel, rmsT0);
rocRMS->SetValue(iChannel, cogRMS);
rocOut->SetValue(iChannel, cogOut);
if ( GetStreamLevel() > 0 ){
TTreeSRedirector *streamer=GetDebugStreamer();
if ( streamer ) {
while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
pad = iChannel-fROC->GetRowIndexes(iSec)[row];
padc = pad-(fROC->GetNPads(iSec,row)/2);
(*streamer) << "DataEnd" <<
"Sector=" << iSec <<
"Pad=" << pad <<
"PadC=" << padc <<
"Row=" << row <<
"PadSec=" << iChannel <<
"Q=" << cogQ <<
"T0=" << cogTime0 <<
"RMS=" << cogRMS <<
"\n";
}
}
}
}
if ( channelCounter>0 ){
fMeanT0rms/=channelCounter;
fMeanQrms/=channelCounter;
fMeanRMSrms/=channelCounter;
}
fVEventTime.ResizeTo(fNevents);
fVEventNumber.ResizeTo(fNevents);
fVTime0SideA.ResizeTo(fNevents);
fVTime0SideC.ResizeTo(fNevents);
}
if (fProcessNew&&fAnalyseNew){
AnalyseTrack();
for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
}
}
}
void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
{
fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
FindLaserLayers();
THnSparse *hProj=fHnDrift;
Double_t posCE[4]={0.,0.,0.,0.};
Double_t widthCE[4]={0.,0.,0.,0.};
for (Int_t i=0; i<4; ++i){
Int_t ce=(i/2>0)*7+6;
hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
TH1 *h=fHnDrift->Projection(3);
h->GetXaxis()->SetRangeUser(fPeaks[ce]-fPeakWidths[ce],fPeaks[ce]+fPeakWidths[ce]);
Int_t nbinMax=h->GetMaximumBin();
Double_t maximum=h->GetMaximum();
Double_t xbinMax=h->GetBinCenter(nbinMax);
TF1 fgaus("gaus","gaus",xbinMax-5,xbinMax+5);
fgaus.SetParameters(maximum,xbinMax,2);
fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
fgaus.SetParLimits(2,0.2,4.);
h->Fit(&fgaus,"RQN");
delete h;
posCE[i]=fgaus.GetParameter(1);
widthCE[i]=4*fgaus.GetParameter(2);
hProj->GetAxis(0)->SetRangeUser(0,72);
}
Float_t vdrift=2.61301900000000000e+06;
AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
Int_t coord[5];
for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
Double_t adc=hProj->GetBinContent(ichunk,coord);
Int_t sector = coord[0]-1;
Int_t row = coord[1]-1;
Int_t pad = coord[2]-1;
Int_t timeBin= coord[3]-1;
Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
Int_t side = (sector/18)%2;
if (time<timestamp-120||time>timestamp+120) continue;
if (adc < 5 ) continue;
if (IsEdgePad(sector,row,pad)) continue;
Int_t padmin=-2, padmax=2;
Int_t timemin=-2, timemax=2;
Int_t minsumperpad=2;
Bool_t isCE=kFALSE;
if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
isCE=kTRUE;
padmin=0;
padmax=0;
timemin=-3;
timemax=7;
}
Bool_t isMaximum=kTRUE;
Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
Double_t cogY=0, rmsY=0;
Int_t npart=0;
for(Int_t ipad=padmin;ipad<=padmax;++ipad){
Float_t lxyz[3];
fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
for(Int_t itime=timemin;itime<=timemax;++itime){
Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
Double_t val=hProj->GetBinContent(a);
totalmass+=val;
tbcm +=(timeBin+itime)*val;
padcm+=(pad+ipad)*val;
cogY +=lxyz[1]*val;
rmstb +=(timeBin+itime)*(timeBin+itime)*val;
rmspad+=(pad+ipad)*(pad+ipad)*val;
rmsY +=lxyz[1]*lxyz[1]*val;
if (val>0) ++npart;
if (val>adc) {
isMaximum=kFALSE;
break;
}
}
if (!isMaximum) break;
}
if (!isMaximum||!npart) continue;
if (totalmass<npart*minsumperpad) continue;
if (!isCE&&rmspad<.1) continue;
tbcm/=totalmass;
padcm/=totalmass;
cogY/=totalmass;
rmstb/=totalmass;
rmspad/=totalmass;
rmsY/=totalmass;
rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
Int_t cog=TMath::Nint(padcm);
Float_t zlength=fParam->GetZLength(side);
Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
Float_t padlxyz[3];
fROC->GetPositionLocal(sector,row,pad,padlxyz);
Double_t gxyz[3]={padlxyz[0],cogY,gz};
Double_t lxyz[3]={padlxyz[0],cogY,gz};
Double_t igxyz[3]={0,0,0};
AliTPCTransform t1;
t1.RotatedGlobal2Global(sector,gxyz);
Double_t mindist=0;
Int_t trackID=-1;
Int_t trackID2=-1;
Int_t index=0;
if (!isCE){
index=row+(sector>35)*fROC->GetNRows(0);
trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
} else {
trackID=336+((sector/18)%2);
index= fROC->GetRowIndexes(sector)[row]+pad;
if (sector<36) {
index+=(sector%18)*fROC->GetNChannels(sector);
} else {
index+=18*fROC->GetNChannels(0);
index+=(sector%18)*fROC->GetNChannels(sector);
}
mindist=1.;
}
if (trackID>0){
AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
Double_t raylength=ltr->GetRayLength();
Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
ltr->GetXYZ(globmir);
if(trackID<336){
if(side==0){
gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
+(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
+(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
}
else {
gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
+(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
+(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
}
}
if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
ltr->fVecSec->GetMatrixArray()[index]=sector;
ltr->fVecP2->GetMatrixArray()[index]=0;
ltr->fVecPhi->GetMatrixArray()[index]=mindist;
ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
}
TObjArray *arr=AliTPCLaserTrack::GetTracks();
ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
}
if (fStreamLevel>4){
(*GetDebugStreamer()) << "clusters" <<
"run=" << fRunNumber <<
"timestamp=" << timestamp <<
"burst=" << burst <<
"side=" << side <<
"sec=" << sector <<
"row=" << row <<
"pad=" << pad <<
"padCog=" << cog <<
"timebin=" << timeBin <<
"cogCE=" << posCE[sector/18] <<
"withCE=" << widthCE[sector/18] <<
"index=" << index <<
"padcm=" << padcm <<
"rmspad=" << rmspad <<
"cogtb=" << tbcm <<
"rmstb=" << rmstb <<
"npad=" << npart <<
"lx=" << padlxyz[0]<<
"ly=" << cogY <<
"lypad=" << padlxyz[1]<<
"rmsY=" << rmsY <<
"gx=" << gxyz[0] <<
"gy=" << gxyz[1] <<
"gz=" << gxyz[2] <<
"igx=" << igxyz[0] <<
"igy=" << igxyz[1] <<
"igz=" << igxyz[2] <<
"mind=" << mindist <<
"max=" << adc <<
"trackid=" << trackID <<
"trackid2=" << trackID2 <<
"npart=" << npart <<
"\n";
}
}
}
void AliTPCCalibCE::AnalyseTrack()
{
AliTPCLaserTrack::LoadTracks();
TObjArray* arrMeasured = SetupMeasured();
TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
AddCEtoIdeal(arrIdeal);
for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
Double_t timestamp=fTimeBursts[iburst];
AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
if (!fHnDrift) continue;
UInt_t entries=(UInt_t)fHnDrift->GetEntries();
if (fBinsLastAna[iburst]>=entries) continue;
fBinsLastAna[iburst]=entries;
for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
ResetMeasured(arrMeasured);
FindLocalMaxima(arrMeasured, timestamp, iburst);
CalculateDV(arrIdeal,arrMeasured,iburst);
if (fStreamLevel>2){
for (Int_t itrack=0; itrack<338; ++itrack){
TObject *iltr=arrIdeal->UncheckedAt(itrack);
TObject *mltr=arrMeasured->UncheckedAt(itrack);
(*GetDebugStreamer()) << "tracks" <<
"run=" << fRunNumber <<
"time=" << timestamp <<
"burst="<< iburst <<
"iltr.=" << iltr <<
"mltr.=" << mltr <<
"\n";
}
}
}
if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
}
Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
const Double_t *peakposloc, Int_t &itrackMin2)
{
TObjArray *arr=AliTPCLaserTrack::GetTracks();
TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
TVector3 vDir;
TVector3 vSt;
Int_t firstbeam=0;
Int_t lastbeam=336/2;
if ( (sector/18)%2 ) {
firstbeam=336/2;
lastbeam=336;
}
mindist=1000000;
Int_t itrackMin=-1;
for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);
vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
Double_t deltaZ=ltr->GetZ()-peakpos[2];
if (TMath::Abs(deltaZ)>40) continue;
vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
vSt.RotateZ(ltr->GetAlpha());
vDir.RotateZ(ltr->GetAlpha());
Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
if (dist<mindist){
mindist=dist;
itrackMin=itrack;
}
}
itrackMin2=-1;
Float_t mindist2=10;
for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack);
if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
Double_t deltaZ=ltr->GetZ()-peakpos[2];
if (TMath::Abs(deltaZ)>40) continue;
Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
if (dist>1) continue;
if (dist<mindist2){
mindist2=dist;
itrackMin2=itrack;
}
}
mindist=mindist2;
return itrackMin2;
}
Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
{
Int_t edge1 = 0;
if ( pad == edge1 ) return kTRUE;
Int_t edge2 = fROC->GetNPads(sector,row)-1;
if ( pad == edge2 ) return kTRUE;
return kFALSE;
}
TObjArray* AliTPCCalibCE::SetupMeasured()
{
TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
TObjArray *arrMeasured = new TObjArray(338);
arrMeasured->SetOwner();
for(Int_t itrack=0;itrack<336;itrack++){
arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
}
AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
ltrce->SetId(336);
ltrce->SetSide(0);
ltrce->fVecSec=new TVectorD(557568/2);
ltrce->fVecP2=new TVectorD(557568/2);
ltrce->fVecPhi=new TVectorD(557568/2);
ltrce->fVecGX=new TVectorD(557568/2);
ltrce->fVecGY=new TVectorD(557568/2);
ltrce->fVecGZ=new TVectorD(557568/2);
ltrce->fVecLX=new TVectorD(557568/2);
ltrce->fVecLY=new TVectorD(557568/2);
ltrce->fVecLZ=new TVectorD(557568/2);
arrMeasured->AddAt(ltrce,336);
ltrce=new AliTPCLaserTrack;
ltrce->SetId(337);
ltrce->SetSide(1);
ltrce->fVecSec=new TVectorD(557568/2);
ltrce->fVecP2=new TVectorD(557568/2);
ltrce->fVecPhi=new TVectorD(557568/2);
ltrce->fVecGX=new TVectorD(557568/2);
ltrce->fVecGY=new TVectorD(557568/2);
ltrce->fVecGZ=new TVectorD(557568/2);
ltrce->fVecLX=new TVectorD(557568/2);
ltrce->fVecLY=new TVectorD(557568/2);
ltrce->fVecLZ=new TVectorD(557568/2);
arrMeasured->AddAt(ltrce,337);
return arrMeasured;
}
void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
{
for(Int_t itrack=0;itrack<338;itrack++){
AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
ltr->fVecSec->Zero();
ltr->fVecP2->Zero();
ltr->fVecPhi->Zero();
ltr->fVecGX->Zero();
ltr->fVecGY->Zero();
ltr->fVecGZ->Zero();
ltr->fVecLX->Zero();
ltr->fVecLY->Zero();
ltr->fVecLZ->Zero();
}
}
void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
{
arr->Expand(338);
AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
ltrceA->SetId(336);
ltrceA->SetSide(0);
ltrceA->fVecSec=new TVectorD(557568/2);
ltrceA->fVecP2=new TVectorD(557568/2);
ltrceA->fVecPhi=new TVectorD(557568/2);
ltrceA->fVecGX=new TVectorD(557568/2);
ltrceA->fVecGY=new TVectorD(557568/2);
ltrceA->fVecGZ=new TVectorD(557568/2);
ltrceA->fVecLX=new TVectorD(557568/2);
ltrceA->fVecLY=new TVectorD(557568/2);
ltrceA->fVecLZ=new TVectorD(557568/2);
arr->AddAt(ltrceA,336);
AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
ltrceC->SetId(337);
ltrceC->SetSide(1);
ltrceC->fVecSec=new TVectorD(557568/2);
ltrceC->fVecP2=new TVectorD(557568/2);
ltrceC->fVecPhi=new TVectorD(557568/2);
ltrceC->fVecGX=new TVectorD(557568/2);
ltrceC->fVecGY=new TVectorD(557568/2);
ltrceC->fVecGZ=new TVectorD(557568/2);
ltrceC->fVecLX=new TVectorD(557568/2);
ltrceC->fVecLY=new TVectorD(557568/2);
ltrceC->fVecLZ=new TVectorD(557568/2);
arr->AddAt(ltrceC,337);
Float_t gxyz[3];
Float_t lxyz[3];
Int_t channelSideA=0;
Int_t channelSideC=0;
Int_t channelSide=0;
AliTPCLaserTrack *ltrce=0x0;
for (Int_t isector=0; isector<72; ++isector){
Int_t side=((isector/18)%2);
for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
fROC->GetPositionLocal(isector,irow,ipad,lxyz);
if (side==0) {
ltrce=ltrceA;
channelSide=channelSideA;
} else {
ltrce=ltrceC;
channelSide=channelSideC;
}
ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
if (side==0){
ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
++channelSideA;
}
else {
ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
++channelSideC;
}
}
}
}
}
void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
{
if (!fArrFitGraphs){
fArrFitGraphs=new TObjArray;
}
Double_t timestamp=fTimeBursts[burst];
static TLinearFitter fdriftA(4,"hyp3");
static TLinearFitter fdriftC(4,"hyp3");
static TLinearFitter fdriftAC(5,"hyp4");
TVectorD fitA(7),fitC(7),fitAC(8);
Float_t chi2A = 10;
Float_t chi2C = 10;
Float_t chi2AC = 10;
Int_t npointsA=0;
Int_t npointsC=0;
Int_t npointsAC=0;
Double_t minres[3]={20.,1,0.8};
for(Int_t i=0;i<3;i++){
fdriftA.ClearPoints();
fdriftC.ClearPoints();
fdriftAC.ClearPoints();
chi2A = 10;
chi2C = 10;
chi2AC = 10;
npointsA=0;
npointsC=0;
npointsAC=0;
for (Int_t itrack=0; itrack<338; ++itrack){
AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
Int_t clustercounter=0;
Int_t indexMax=159;
for (Int_t index=0; index<indexMax; ++index){
Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
}
if (clustercounter>130&&itrack<336) continue;
clustercounter=0;
Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
if (itrack>335) indexMax=557568/2;
for (Int_t index=0; index<indexMax; ++index){
Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
if (iltr->GetBundle()==0) continue;
if (iR<133||mR<133) continue;
if(TMath::Abs(mltr->fVecP2->GetMatrixArray()[index])>minres[i]) continue;
Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
if (iltr->GetSide()==0){
fdriftA.AddPoint(xxx,mdrift,1);
}else{
fdriftC.AddPoint(xxx,mdrift,1);
}
Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, static_cast<Double_t>(iltr->GetSide())};
fdriftAC.AddPoint(xxx2,mdrift,1);
}
}
fdriftA.Eval();
fdriftC.Eval();
fdriftAC.Eval();
fdriftA.GetParameters(fitA);
fdriftC.GetParameters(fitC);
fdriftAC.GetParameters(fitAC);
for (Int_t itrack=0; itrack<338; ++itrack){
AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
Int_t clustercounter=0;
Int_t indexMax=159;
for (Int_t index=0; index<indexMax; ++index){
Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
}
if (clustercounter>130&&itrack<336) continue;
clustercounter=0;
Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
if (itrack>335) indexMax=557568/2;
for (Int_t index=0; index<indexMax; ++index){
Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
if (iR<60||mR<60) continue;
Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
TVectorD *v=&fitA;
if (iltr->GetSide()==1) v=&fitC;
Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
}
}
fitA.ResizeTo(7);
fitC.ResizeTo(7);
fitAC.ResizeTo(8);
npointsA= fdriftA.GetNpoints();
if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
fitA[5]=npointsA;
fitA[6]=chi2A;
npointsC= fdriftC.GetNpoints();
if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
fitC[5]=npointsC;
fitC[6]=chi2C;
npointsAC= fdriftAC.GetNpoints();
if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
fitAC[5]=npointsAC;
fitAC[6]=chi2AC;
if (fStreamLevel>2){
(*GetDebugStreamer()) << "DriftV" <<
"iter=" << i <<
"run=" << fRunNumber <<
"time=" << timestamp <<
"fitA.=" << &fitA <<
"fitC.=" << &fitC <<
"fitAC.=" << &fitAC <<
"\n";
}
}
TGraphErrors *grA[7];
TGraphErrors *grC[7];
TGraphErrors *grAC[8];
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_");
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_");
TObjArray *arrNames=names.Tokenize(";");
TObjArray *arrNamesAC=namesAC.Tokenize(";");
for (Int_t i=0; i<7; ++i){
TString grName=arrNames->UncheckedAt(i)->GetName();
grName+="A";
grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
if (!grA[i]){
grA[i]=new TGraphErrors;
grA[i]->SetName(grName.Data());
grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
fArrFitGraphs->Add(grA[i]);
}
Int_t ipoint=burst;
grA[i]->SetPoint(ipoint,timestamp,fitA(i));
grA[i]->SetPointError(ipoint,60,.0001);
if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
}
for (Int_t i=0; i<7; ++i){
TString grName=arrNames->UncheckedAt(i)->GetName();
grName+="C";
grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
if (!grC[i]){
grC[i]=new TGraphErrors;
grC[i]->SetName(grName.Data());
grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
fArrFitGraphs->Add(grC[i]);
}
Int_t ipoint=burst;
grC[i]->SetPoint(ipoint,timestamp,fitC(i));
grC[i]->SetPointError(ipoint,60,.0001);
if (i<4) grC[i]->SetPointError(ipoint,60,fdriftC.GetCovarianceMatrixElement(i,i));
}
for (Int_t i=0; i<8; ++i){
TString grName=arrNamesAC->UncheckedAt(i)->GetName();
grName+="AC";
grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
if (!grAC[i]){
grAC[i]=new TGraphErrors;
grAC[i]->SetName(grName.Data());
grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
fArrFitGraphs->Add(grAC[i]);
}
Int_t ipoint=burst;
grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
grAC[i]->SetPointError(ipoint,60,.0001);
if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftAC.GetCovarianceMatrixElement(i,i));
}
if (fDebugLevel>10){
printf("A side fit parameters:\n");
fitA.Print();
printf("\nC side fit parameters:\n");
fitC.Print();
printf("\nAC side fit parameters:\n");
fitAC.Print();
}
delete arrNames;
delete arrNamesAC;
}
Double_t AliTPCCalibCE::SetBurstHnDrift()
{
Int_t i=0;
for(i=0; i<fTimeBursts.GetNrows(); ++i){
if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
return fTimeBursts(i);
}
}
CreateDVhist();
fArrHnDrift.AddAt(fHnDrift,i);
fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
for (i=0;i<14;++i){
fPeaks[i]=0;
fPeakWidths[i]=0;
}
fEventInBunch=0;
return fTimeStamp;
}
void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t )
{
TString sDir(dir);
TString objName=GetName();
Int_t type=0;
TObjArray *arr=sDir.Tokenize(",");
TIter next(arr);
TObjString *s;
while ( (s=(TObjString*)next()) ){
TString optString=s->GetString();
optString.Remove(TString::kBoth,' ');
if (optString.BeginsWith("name=")){
objName=optString.ReplaceAll("name=","");
}
if (optString.BeginsWith("type=")){
optString.ReplaceAll("type=","");
type=optString.Atoi();
}
}
delete arr;
if ( type==4 ){
AliTPCCalibCE ce;
ce.fArrFitGraphs=fArrFitGraphs;
ce.fNevents=fNevents;
ce.fTimeBursts.ResizeTo(fTimeBursts.GetNrows());
ce.fTimeBursts=fTimeBursts;
ce.fProcessNew=kTRUE;
TFile f(filename,"recreate");
ce.Write(objName.Data());
fArrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
f.Close();
ce.fArrFitGraphs=0x0;
return;
}
if (type==1||type==2) {
fCalRocArrayT0.Delete();
fCalRocArrayT0Err.Delete();
fCalRocArrayQ.Delete();
fCalRocArrayRMS.Delete();
fCalRocArrayOutliers.Delete();
}
if (type==2||type==3){
fParamArrayEventPol1.Delete();
fParamArrayEventPol2.Delete();
}
TObjArray histoQArray(72);
TObjArray histoT0Array(72);
TObjArray histoRMSArray(72);
TObjArray arrHnDrift(fArrHnDrift.GetEntries());
if (type==1||type==2||type==3){
for (Int_t i=0; i<72; ++i){
histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
}
fHistoQArray.SetOwner(kFALSE);
fHistoT0Array.SetOwner(kFALSE);
fHistoRMSArray.SetOwner(kFALSE);
fHistoQArray.Clear();
fHistoT0Array.Clear();
fHistoRMSArray.Clear();
for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
}
fArrHnDrift.SetOwner(kFALSE);
fArrHnDrift.Clear();
}
TDirectory *backup = gDirectory;
TFile f(filename,"recreate");
Write(objName.Data());
if (type==1||type==2) {
histoQArray.Write("histoQArray",TObject::kSingleKey);
histoT0Array.Write("histoT0Array",TObject::kSingleKey);
histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
}
f.Save();
f.Close();
if (type==1||type==2){
for (Int_t i=0; i<72; ++i){
fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
}
fHistoQArray.SetOwner(kTRUE);
fHistoT0Array.SetOwner(kTRUE);
fHistoRMSArray.SetOwner(kTRUE);
for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
}
fArrHnDrift.SetOwner(kTRUE);
}
if ( backup ) backup->cd();
}
AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
{
TFile f(filename);
if (!f.IsOpen() || f.IsZombie() ) return 0x0;
TList *l=f.GetListOfKeys();
TIter next(l);
TKey *key=0x0;
TObject *o=0x0;
AliTPCCalibCE *ce=0x0;
TObjArray *histoQArray=0x0;
TObjArray *histoT0Array=0x0;
TObjArray *histoRMSArray=0x0;
TObjArray *arrHnDrift=0x0;
while ( (key=(TKey*)next()) ){
o=key->ReadObj();
if ( o->IsA()==AliTPCCalibCE::Class() ){
ce=(AliTPCCalibCE*)o;
} else if ( o->IsA()==TObjArray::Class() ){
TString name=key->GetName();
if ( name=="histoQArray") histoQArray=(TObjArray*)o;
if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
}
}
if (ce){
TH1* hist=0x0;
if (histoQArray){
for (Int_t i=0; i<72; ++i){
hist=(TH1*)histoQArray->UncheckedAt(i);
if (hist) hist->SetDirectory(0x0);
ce->fHistoQArray.AddAt(hist,i);
}
ce->fHistoQArray.SetOwner(kTRUE);
}
if (histoT0Array) {
for (Int_t i=0; i<72; ++i){
hist=(TH1*)histoT0Array->UncheckedAt(i);
if (hist) hist->SetDirectory(0x0);
ce->fHistoT0Array.AddAt(hist,i);
}
ce->fHistoT0Array.SetOwner(kTRUE);
}
if (histoRMSArray){
for (Int_t i=0; i<72; ++i){
hist=(TH1*)histoRMSArray->UncheckedAt(i);
if (hist) hist->SetDirectory(0x0);
ce->fHistoRMSArray.AddAt(hist,i);
}
ce->fHistoRMSArray.SetOwner(kTRUE);
}
if (arrHnDrift){
for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
ce->fArrHnDrift.AddAt(hSparse,i);
}
}
ce->Analyse();
}
f.Close();
return ce;
}