From d0b31f577d4524c9fbbd8d88e76c4b1df3ba743b Mon Sep 17 00:00:00 2001 From: marian Date: Tue, 28 Apr 2009 15:27:32 +0000 Subject: [PATCH] Bug removal Adding additional functionality - example to create the spline fit --- TPC/AliTPCcalibBase.cxx | 12 ++- TPC/AliTPCcalibTimeGain.cxx | 198 +++++++++++++++++++----------------- TPC/AliTPCcalibTimeGain.h | 22 ++-- 3 files changed, 128 insertions(+), 104 deletions(-) diff --git a/TPC/AliTPCcalibBase.cxx b/TPC/AliTPCcalibBase.cxx index fa33ab7331e..b6f6e0d9216 100644 --- a/TPC/AliTPCcalibBase.cxx +++ b/TPC/AliTPCcalibBase.cxx @@ -283,15 +283,19 @@ TGraphErrors * AliTPCcalibBase::FitSlices(THnSparse *h, Int_t axisDim1, Int_t ax Float_t xMin = projectionHist->GetXaxis()->GetXmin(); Float_t xMax = projectionHist->GetXaxis()->GetXmax(); Float_t xMedian = (xMin+xMax)*0.5; - Float_t integral = projectionHist->GetSum(); + Float_t integral = 0; + for(Int_t jbin=1; jbinGetNbinsX()-1; jbin++) { + integral+=projectionHist->GetBinContent(jbin); + } + printf("Integral %f\t%f\n",integral, projectionHist->GetSum()); // // Float_t currentSum=0; for(Int_t jbin=1; jbinGetNbinsX()-1; jbin++) { currentSum += projectionHist->GetBinContent(jbin); - if (currentSum/integralGetBinCenter(jbin); - if (currentSum/integralGetBinCenter(jbin+1); - if (currentSum/integral<0.5 && projectionHist->GetBinContent(jbin)>0){ + if (currentSumGetBinCenter(jbin); + if (currentSumGetBinCenter(jbin+1); + if (currentSum<0.5*integral && projectionHist->GetBinContent(jbin)>0){ xMedian = (projectionHist->GetBinCenter(jbin)*projectionHist->GetBinContent(jbin)+ projectionHist->GetBinCenter(jbin+1)*projectionHist->GetBinContent(jbin+1))/ (projectionHist->GetBinContent(jbin)+projectionHist->GetBinContent(jbin+1)); diff --git a/TPC/AliTPCcalibTimeGain.cxx b/TPC/AliTPCcalibTimeGain.cxx index cc0fe79d103..96c9669ae99 100644 --- a/TPC/AliTPCcalibTimeGain.cxx +++ b/TPC/AliTPCcalibTimeGain.cxx @@ -26,13 +26,15 @@ gSystem->Load("libTPCcalib"); // compare reference // -//2. Visulaize results +//2. Visualize results // TFile fcalib("CalibObjects.root"); TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); AliTPCcalibTimeGain * gain = ( AliTPCcalibTimeGain *)array->FindObject("calibTimeGain"); -TGraph * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10); -gr->Draw("ALPsame") +TGraphErrors * gr = AliTPCcalibBase::FitSlices(gain->GetHistGainTime(),0,1,2000,10,0.2,0.8); +gain->GetHistGainTime()->Projection(0,1)->Draw("colz") +gr->SetMarkerStyle(25); +gr->Draw("lp") // @@ -61,6 +63,7 @@ grfit->Draw("lu"); #include "TH2F.h" #include "TH3F.h" #include "THnSparse.h" +#include "TGraphErrors.h" #include "TList.h" #include "TMath.h" #include "TCanvas.h" @@ -110,7 +113,8 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain() fUseShapeNorm(0), fUsePosNorm(0), fUsePadNorm(0), - fIsCosmic(0) + fIsCosmic(0), + fLowMemoryConsumption(0) { AliInfo("Default Constructor"); } @@ -128,7 +132,8 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title fUseShapeNorm(0), fUsePosNorm(0), fUsePadNorm(0), - fIsCosmic(0) + fIsCosmic(0), + fLowMemoryConsumption(0) { SetName(name); @@ -141,11 +146,12 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title Double_t deltaTime = EndTime - StartTime; - // main histogram for time dependence: dE/dx, time, type (0-all,1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available) - Int_t binsGainTime[5] = {100, TMath::Nint(deltaTime/deltaIntegrationTimeGain), 3, 25, 200}; - Double_t xminGainTime[5] = {0.5, StartTime, -0.5, 0, 0.1}; + // main histogram for time dependence: dE/dx, time, type (1-muon cosmic,2-pion beam data), meanDriftlength, momenta (only filled if enough space is available) + Int_t timeBins = TMath::Nint(deltaTime/deltaIntegrationTimeGain); + Int_t binsGainTime[5] = {100, timeBins, 2, 25, 200}; + Double_t xminGainTime[5] = {0.5, StartTime, 0.5, 0, 0.1}; Double_t xmaxGainTime[5] = { 4, EndTime, 2.5, 250, 50}; - fHistGainTime = new THnSparseF("HistGainTime","dEdx l;time;dEdx",5,binsGainTime,xminGainTime,xmaxGainTime); + fHistGainTime = new THnSparseF("HistGainTime","dEdx time dep.;dEdx;dEdx,time,type,driftlength,momenta",5,binsGainTime,xminGainTime,xmaxGainTime); BinLogX(fHistGainTime, 4); // fHistDeDxTotal = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",250,0.01,100.,1000,0.,1000); @@ -160,6 +166,7 @@ AliTPCcalibTimeGain::AliTPCcalibTimeGain(const Text_t *name, const Text_t *title fUsePadNorm = kFALSE; // fIsCosmic = kTRUE; + fLowMemoryConsumption = kTRUE; // } @@ -180,10 +187,22 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) { if (!event) { Printf("ERROR: ESD not available"); return; - } + } + + if (fIsCosmic) { // this should be removed at some point based on trigger mask !? + ProcessCosmicEvent(event); + } else { + ProcessBeamEvent(event); + } + - Int_t ntracks=event->GetNumberOfTracks(); + +} + + +void AliTPCcalibTimeGain::ProcessCosmicEvent(AliESDEvent *event) { + AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); if (!ESDfriend) { Printf("ERROR: ESDfriend not available"); @@ -191,7 +210,7 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) { } // UInt_t time = event->GetTimeStamp(); - if (time < 0.1) time = (UInt_t)(fHistGainTime->GetAxis(1)->GetXmin() + 1); + Int_t ntracks = event->GetNumberOfTracks(); // // track loop // @@ -206,23 +225,12 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) { if (!trackOut) continue; // calculate necessary track parameters - //Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP()); Double_t meanP = trackIn->GetP(); Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ()); - //Double_t d = trackIn->GetLinearD(0,0); Int_t NclsDeDx = track->GetTPCNcls(); - //if (meanP > 0.7 || meanP < 0.2) continue; - if (fIsCosmic && meanP < 20) continue; - if (NclsDeDx < 60) continue; - // exclude tracks which do not look like primaries or are simply too short or on wrong sectors - - //if (TMath::Abs(trackIn->GetSnp()) > 3*0.4) continue; - //if (TMath::Abs(trackIn->GetZ()) > 150) continue; - //if (seed->CookShape(1) > 1) continue; - //if (TMath::Abs(trackIn->GetY()) > 20) continue; - //if (TMath::Abs(d)>20) continue; // distance to the 0,0; select only tracks which cross chambers under proper angle + if (NclsDeDx < 60) continue; if (TMath::Abs(trackIn->GetTgl()) > 1) continue; if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue; @@ -237,31 +245,80 @@ void AliTPCcalibTimeGain::Process(AliESDEvent *event) { if (seed) { Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0); - //Double_t TPCsignalMax = (1/fMIP)*track->GetTPCsignal(); fHistDeDxTotal->Fill(meanP, TPCsignalMax); - // - //dE/dx, time, type (0-all,1-muon cosmic,2-pion beam data), momenta - Double_t vec[5] = {TPCsignalMax,time,0,meanDrift,meanP}; - fHistGainTime->Fill(vec); // avoid this filling if memory consumption is too high - - // only partial filling if memory consumption has to be kept low; for cosmic and beam data - if (fIsCosmic) { - Double_t vecCos[5] = {TPCsignalMax,time,1,meanDrift,20}; - if (meanP > 20) fHistGainTime->Fill(vecCos); - } else { - Double_t vecBeam[5] = {TPCsignalMax,time,2,meanDrift,0.5}; - if (meanP > 0.4 && meanP < 0.5) fHistGainTime->Fill(vecBeam); + if (fLowMemoryConsumption) { + if (meanP < 20) continue; + meanP = 30; // set all momenta to one in order to save memory } + //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta + Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP}; + fHistGainTime->Fill(vec); + } + } - } else { - cout << "ERROR: TPC seed not found" << endl; +} + + + +void AliTPCcalibTimeGain::ProcessBeamEvent(AliESDEvent *event) { + + AliESDfriend *ESDfriend=static_cast(event->FindListObject("AliESDfriend")); + if (!ESDfriend) { + Printf("ERROR: ESDfriend not available"); + return; + } + // + UInt_t time = event->GetTimeStamp(); + Int_t ntracks = event->GetNumberOfTracks(); + // + // track loop + // + for (Int_t i=0;iGetTrack(i); + if (!track) continue; + + const AliExternalTrackParam * trackIn = track->GetInnerParam(); + const AliExternalTrackParam * trackOut = track->GetOuterParam(); + if (!trackIn) continue; + if (!trackOut) continue; + + // calculate necessary track parameters + Double_t meanP = trackIn->GetP(); + Double_t meanDrift = 250 - 0.5*TMath::Abs(trackIn->GetZ() + trackOut->GetZ()); + Int_t NclsDeDx = track->GetTPCNcls(); + + // exclude tracks which do not look like primaries or are simply too short or on wrong sectors + if (NclsDeDx < 60) continue; + if (TMath::Abs(trackIn->GetTgl()) > 1) continue; + if (TMath::Abs(trackIn->GetSnp()) > 0.6) continue; + + // Get seeds + AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i); + if (!friendTrack) continue; + TObject *calibObject; + AliTPCseed *seed = 0; + for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) { + if ((seed=dynamic_cast(calibObject))) break; + } + + if (seed) { + Double_t TPCsignalMax = (1/fMIP)*seed->CookdEdxNorm(fLowerTrunc,fUpperTrunc,1,0,159,fUseShapeNorm,fUsePosNorm,fUsePadNorm,0); + fHistDeDxTotal->Fill(meanP, TPCsignalMax); + // + if (fLowMemoryConsumption) { + if (meanP > 0.5 || meanP < 0.4) continue; + meanP = 0.45; // set all momenta to one in order to save memory + } + //dE/dx, time, type (1-muon cosmic,2-pion beam data), momenta + Double_t vec[5] = {TPCsignalMax,time,1,meanDrift,meanP}; + fHistGainTime->Fill(vec); } } - - + } @@ -269,16 +326,15 @@ void AliTPCcalibTimeGain::Analyze() { // // // - TObjArray arr; if (fIsCosmic) { - fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); + fHistGainTime->GetAxis(2)->SetRangeUser(0.51,1.49); // only cosmics + fHistGainTime->GetAxis(4)->SetRangeUser(20,100); // only Fermi-Plateau muons } else { - fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); + fHistGainTime->GetAxis(2)->SetRangeUser(1.51,2.49); // only beam data + fHistGainTime->GetAxis(4)->SetRangeUser(0.39,0.51); // only MIP pions } - fHistGainTime->Projection(0,1)->FitSlicesY(0,0,-1,0,"QNR",&arr); - TH1D * fitMean = (TH1D*) arr.At(1); // - fGainVsTime = new TGraph(fitMean); + fGainVsTime = AliTPCcalibBase::FitSlices(fHistGainTime,0,1,2000,10); // return; } @@ -306,48 +362,6 @@ Long64_t AliTPCcalibTimeGain::Merge(TCollection *li) { } -TGraph * AliTPCcalibTimeGain::FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries){ - // - // Fitting slices of the projection(axisDim1,axisDim2) of a sparse histogram - // - TH2D * hist = h->Projection(axisDim1, axisDim2); - - Double_t xvec[10000]; - Double_t yvec[10000]; - Int_t counter = 0; - - for(Int_t i=1; i < hist->GetNbinsX(); i++) { - Int_t interval = 0; - if (hist->Integral(i,i,0,hist->GetNbinsY()) < minEntries) { - if (hist->Integral(i,i+1,0,hist->GetNbinsY()) < minEntries) { - if (hist->Integral(i,i+2,0,hist->GetNbinsY()) < minEntries) { - continue; - } else { - interval = 2; - } - } else { - interval = 1; - } - } - counter++; - i += interval; - // - Double_t x = hist->GetXaxis()->GetBinCenter(i); - TH1D * projectionHist = hist->ProjectionY("dummy",i,i + interval); - TF1 funcGaus("funcGaus","gaus"); - projectionHist->Fit(&funcGaus,"QN"); - // - xvec[counter-1] = x; - yvec[counter-1] = funcGaus.GetParameter(1); - delete projectionHist; - } - - TGraph * graph = new TGraph(counter, xvec, yvec); - - return graph; -} - - void AliTPCcalibTimeGain::BinLogX(THnSparse *h, Int_t axisDim) { @@ -400,10 +414,10 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) //{0.0762*MIP,10.632,1.34e-05,1.863,1.948} const Double_t sigma = 0.06; - TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",hist->GetNbinsX(),0.1,5000.,hist->GetNbinsY(),0.5,5.); + TH2F *histBG = new TH2F("histBG","dEdxBg; #beta #gamma; TPC signal (a.u.)",TMath::Nint(0.5*hist->GetNbinsX()),0.1,5000.,hist->GetNbinsY(),hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax()); BinLogX(histBG); - TF1 *foElectron = new TF1("foElectron", "(1.7/1.6)*AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100); + TF1 *foElectron = new TF1("foElectron", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.01,100); TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.1056,[0],[1],[2],[3],[4])",0.01,100); TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.138,[0],[1],[2],[3],[4])",0.01,100); TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.498,[0],[1],[2],[3],[4])",0.01,100); @@ -445,7 +459,7 @@ void AliTPCcalibTimeGain::CalculateBetheAlephParams(TH2F *hist, Double_t * ini) for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y); } } else { - if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 3*sigma && ((x>0.25&&x<0.35) || (x>1.5&&x<1.8) || (x>0.65&&x<0.7))) { + if (TMath::Abs(y - foElectron->Eval(x))/foElectron->Eval(x) < 2*sigma && ((x>0.25&&x<0.35) || (x>1.5&&x<1.8) || (x>0.65&&x<0.7))) { for(Int_t iEntries =0; iEntries < entries; iEntries++) histBG->Fill(x/mass, y); } } diff --git a/TPC/AliTPCcalibTimeGain.h b/TPC/AliTPCcalibTimeGain.h index 5e31facc0e0..c093df4935f 100644 --- a/TPC/AliTPCcalibTimeGain.h +++ b/TPC/AliTPCcalibTimeGain.h @@ -15,13 +15,14 @@ class TH3F; class TH2F; class THnSparse; class TList; +class TGraphErrors; class AliESDEvent; class AliESDtrack; class AliTPCcalibLaser; #include "TTreeStream.h" - + class AliTPCcalibTimeGain:public AliTPCcalibBase { public: AliTPCcalibTimeGain(); @@ -32,15 +33,18 @@ public: virtual Long64_t Merge(TCollection *li); virtual void Analyze(); // + void ProcessCosmicEvent(AliESDEvent *event); + void ProcessBeamEvent(AliESDEvent *event); + // void CalculateBetheAlephParams(TH2F *hist, Double_t * ini); static void BinLogX(THnSparse *h, Int_t axisDim); static void BinLogX(TH1 *h); - static TGraph * FitSlices(THnSparse *h, Int_t axisDim1, Int_t axisDim2, Int_t minEntries); // THnSparse * GetHistGainTime(){return (THnSparse*) fHistGainTime;}; - TGraph * GetGraphGainVsTime(){return fGainVsTime;}; TH2F * GetHistDeDxTotal(){return (TH2F*) fHistDeDxTotal;}; // + TGraphErrors * GetGraphGainVsTime(){return fGainVsTime;}; + // void SetMIP(Float_t MIP){fMIP = MIP;}; void SetLowerTrunc(Float_t LowerTrunc){fLowerTrunc = LowerTrunc;}; void SetUpperTrunc(Float_t UpperTrunc){fUpperTrunc = UpperTrunc;}; @@ -48,16 +52,17 @@ public: void SetUsePosNorm(Bool_t UsePosNorm){fUsePosNorm = UsePosNorm;}; void SetUsePadNorm(Int_t UsePadNorm){fUsePadNorm = UsePadNorm;}; void SetIsCosmic(Bool_t IsCosmic){fIsCosmic = IsCosmic;}; + void SetLowMemoryConsumption(Bool_t LowMemoryConsumption){fLowMemoryConsumption = LowMemoryConsumption;}; private: // - THnSparse * fHistGainTime; // dEdx vs. time, type, Driftlength, momentum P - TGraph * fGainVsTime; // multiplication factor vs. time - TH2F * fHistDeDxTotal; // dEdx vs. momentum for quality assurance + THnSparse * fHistGainTime; // dEdx vs. time, type, Driftlength, momentum P + TGraphErrors * fGainVsTime; // multiplication factor vs. time + TH2F * fHistDeDxTotal; // dEdx vs. momentum for quality assurance // Float_t fIntegrationTimeDeDx; // required statistics for each dEdx time bin // - Float_t fMIP; // rough MIP position in order to have scalable histograms + Float_t fMIP; // rough MIP position in order to have scaleable histograms // Float_t fLowerTrunc; // lower truncation of dE/dx ; at most 5% Float_t fUpperTrunc; // upper truncation of dE/dx ; ca. 70% @@ -65,7 +70,8 @@ private: Bool_t fUsePosNorm; // charge correction (analytical?) Int_t fUsePadNorm; // normalization of pad geometries // - Bool_t fIsCosmic; // kTRUE if the analyzed runs are contain cosmic events + Bool_t fIsCosmic; // kTRUE if the analyzed runs contain cosmic events + Bool_t fLowMemoryConsumption; // set this option kTRUE if the momenta information should not be stored in order to save memory // AliTPCcalibTimeGain(const AliTPCcalibTimeGain&); AliTPCcalibTimeGain& operator=(const AliTPCcalibTimeGain&); -- 2.43.0