X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGLF%2FFORWARD%2Fanalysis2%2FAliCentralMCCorrectionsTask.cxx;h=39c193fff3e04d81db108b3333250ab624af82a4;hb=b8f92f9dec22bca57816e5765e3216b63199ae46;hp=1fe6d1f37d32f2b62f4a19b1bd26b64908e7d73c;hpb=438d5c4ef9c4e8290b7c8195cfda002ee1112dc6;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx b/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx index 1fe6d1f37d3..39c193fff3e 100644 --- a/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx +++ b/PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx @@ -12,6 +12,7 @@ // Corrections used // #include "AliCentralMCCorrectionsTask.h" +#include "AliCentralCorrectionManager.h" #include "AliTriggerAnalysis.h" #include "AliPhysicsSelection.h" #include "AliLog.h" @@ -45,19 +46,14 @@ namespace { //==================================================================== AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask() - : AliAnalysisTaskSE(), - fInspector(), + : AliBaseMCCorrectionsTask(), fTrackDensity(), - fVtxBins(0), - fFirstEvent(true), - fHEvents(0), - fHEventsTr(0), - fHEventsTrVtx(0), - fVtxAxis(), - fEtaAxis(), - fList(), + fSecCorr(0), + fAccCorr(0), fNPhiBins(20), - fEffectiveCorr(true) + fEffectiveCorr(true), + fEtaCut(1.9), + fCorrCut(0.6) { // // Constructor @@ -65,24 +61,19 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask() // Parameters: // name Name of task // - DGUARD(fDebug,0,"Default construction of AliCentralMCCorrectionsTask"); + DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask"); } //____________________________________________________________________ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name) - : AliAnalysisTaskSE(name), - fInspector("eventInspector"), + : AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())), fTrackDensity("trackDensity"), - fVtxBins(0), - fFirstEvent(true), - fHEvents(0), - fHEventsTr(0), - fHEventsTrVtx(0), - fVtxAxis(10,-10,10), - fEtaAxis(200,-4,6), - fList(), + fSecCorr(0), + fAccCorr(0), fNPhiBins(20), - fEffectiveCorr(true) + fEffectiveCorr(true), + fEtaCut(1.9), + fCorrCut(0.6) { // // Constructor @@ -90,432 +81,94 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name) // Parameters: // name Name of task // - DGUARD(fDebug,0,"Named construction of AliCentralMCCorrectionsTask: %s",name); - DefineOutput(1, TList::Class()); - DefineOutput(2, TList::Class()); + DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name); } //____________________________________________________________________ -AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o) - : AliAnalysisTaskSE(o), - fInspector(o.fInspector), - fTrackDensity(), - fVtxBins(0), - fFirstEvent(o.fFirstEvent), - fHEvents(o.fHEvents), - fHEventsTr(o.fHEventsTr), - fHEventsTrVtx(o.fHEventsTrVtx), - fVtxAxis(10,-10,10), - fEtaAxis(200,-4,6), - fList(o.fList), - fNPhiBins(o.fNPhiBins), - fEffectiveCorr(o.fEffectiveCorr) -{ - // - // Copy constructor - // - // Parameters: - // o Object to copy from - // - DGUARD(fDebug,0,"Copy construction of AliCentralMCCorrectionsTask"); -} - -//____________________________________________________________________ -AliCentralMCCorrectionsTask& -AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o) -{ - // - // Assignment operator - // - // Parameters: - // o Object to assign from - // - // Return: - // Reference to this object - // - DGUARD(fDebug,3,"Assignment of AliCentralMCCorrectionsTask"); - if (&o == this) return *this; - fInspector = o.fInspector; - fTrackDensity = o.fTrackDensity; - fVtxBins = o.fVtxBins; - fFirstEvent = o.fFirstEvent; - fHEvents = o.fHEvents; - fHEventsTr = o.fHEventsTr; - fHEventsTrVtx = o.fHEventsTrVtx; - SetVertexAxis(o.fVtxAxis); - SetEtaAxis(o.fEtaAxis); - fNPhiBins = o.fNPhiBins; - fEffectiveCorr = o.fEffectiveCorr; - - return *this; -} - -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::Init() -{ - // - // Initialize the task - // - // - DGUARD(fDebug,1,"Initialize AliCentralMCCorrectionsTask"); -} - -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, - Double_t max) -{ - // - // Set the vertex axis to use - // - // Parameters: - // nBins Number of bins - // vzMin Least @f$z@f$ coordinate of interation point - // vzMax Largest @f$z@f$ coordinate of interation point - // - DGUARD(fDebug,3,"Set vertex axis AliCentralMCCorrectionsTask [%d,%f,%f]", - nBin, min, max); - if (max < min) { - Double_t tmp = min; - min = max; - max = tmp; - } - if (min < -10) - AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min)); - if (max > +10) - AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max)); - fVtxAxis.Set(nBin, min, max); -} -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis) -{ - // - // Set the vertex axis to use - // - // Parameters: - // axis Axis - // - SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax()); -} - -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max) -{ - // - // Set the eta axis to use - // - // Parameters: - // nBins Number of bins - // vzMin Least @f$\eta@f$ - // vzMax Largest @f$\eta@f$ - // - DGUARD(fDebug,3,"Set eta axis AliCentralMCCorrectionsTask [%d,%f,%f]", - nBin, min, max); - if (max < min) { - Double_t tmp = min; - min = max; - max = tmp; - } - if (min < -4) - AliWarning(Form("Minimum eta %f < -4, make sure you want this",min)); - if (max > +6) - AliWarning(Form("Minimum eta %f > +6, make sure you want this",max)); - fEtaAxis.Set(nBin, min, max); -} -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis) +AliBaseMCCorrectionsTask::VtxBin* +AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high) { - // - // Set the eta axis to use - // - // Parameters: - // axis Axis - // - SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax()); + return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins); } //____________________________________________________________________ -void -AliCentralMCCorrectionsTask::DefineBins(TList* l) -{ - DGUARD(fDebug,1,"Define bins in AliCentralMCCorrectionsTask"); - if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1); - if (fVtxBins->GetEntries() > 0) return; - - fVtxBins->SetOwner(); - for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) { - Double_t low = fVtxAxis.GetBinLowEdge(i); - Double_t high = fVtxAxis.GetBinUpEdge(i); - VtxBin* bin = new VtxBin(low, high, fEtaAxis, fNPhiBins); - fVtxBins->AddAt(bin, i); - bin->DefineOutput(l); - } -} - -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::UserCreateOutputObjects() -{ - // - // Create output objects - // - // - DGUARD(fDebug,1,"Create user output for AliCentralMCCorrectionsTask"); - fList = new TList; - fList->SetOwner(); - fList->SetName(Form("%sSums", GetName())); - - DefineBins(fList); - - fHEvents = new TH1I(GetEventName(false,false), - "Number of all events", - fVtxAxis.GetNbins(), - fVtxAxis.GetXmin(), - fVtxAxis.GetXmax()); - fHEvents->SetXTitle("v_{z} [cm]"); - fHEvents->SetYTitle("# of events"); - fHEvents->SetFillColor(kBlue+1); - fHEvents->SetFillStyle(3001); - fHEvents->SetDirectory(0); - fList->Add(fHEvents); - - fHEventsTr = new TH1I(GetEventName(true, false), - "Number of triggered events", - fVtxAxis.GetNbins(), - fVtxAxis.GetXmin(), - fVtxAxis.GetXmax()); - fHEventsTr->SetXTitle("v_{z} [cm]"); - fHEventsTr->SetYTitle("# of events"); - fHEventsTr->SetFillColor(kRed+1); - fHEventsTr->SetFillStyle(3001); - fHEventsTr->SetDirectory(0); - fList->Add(fHEventsTr); - - fHEventsTrVtx = new TH1I(GetEventName(true,true), - "Number of events w/trigger and vertex", - fVtxAxis.GetNbins(), - fVtxAxis.GetXmin(), - fVtxAxis.GetXmax()); - fHEventsTrVtx->SetXTitle("v_{z} [cm]"); - fHEventsTrVtx->SetYTitle("# of events"); - fHEventsTrVtx->SetFillColor(kBlue+1); - fHEventsTrVtx->SetFillStyle(3001); - fHEventsTrVtx->SetDirectory(0); - fList->Add(fHEventsTrVtx); - - // Copy axis objects to output - TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis", - fVtxAxis.GetNbins(), - fVtxAxis.GetXmin(), - fVtxAxis.GetXmax()); - TH1* etaAxis = new TH1D("etaAxis", "Eta axis", - fEtaAxis.GetNbins(), - fEtaAxis.GetXmin(), - fEtaAxis.GetXmax()); - fList->Add(vtxAxis); - fList->Add(etaAxis); - - AliInfo(Form("Initialising sub-routines: %p, %p", - &fInspector, &fTrackDensity)); - fInspector.DefineOutput(fList); - fTrackDensity.DefineOutput(fList); - - PostData(1, fList); -} -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::UserExec(Option_t*) +Bool_t +AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd, + const AliMCEvent& mc, + AliBaseMCCorrectionsTask::VtxBin& bin, + Double_t vz) { - // - // Process each event - // - // Parameters: - // option Not used - // - - DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event"); - // Get the input data - MC event - AliMCEvent* mcEvent = MCEvent(); - if (!mcEvent) { - AliWarning("No MC event found"); - return; - } - - // Get the input data - ESD event - AliESDEvent* esd = dynamic_cast(InputEvent()); - if (!esd) { - AliWarning("No ESD event found for input event"); - return; - } - - //--- Read run information ----------------------------------------- - if (fFirstEvent && esd->GetESDRun()) { - fInspector.ReadRunDetails(esd); - - AliInfo(Form("Initializing with parameters from the ESD:\n" - " AliESDEvent::GetBeamEnergy() ->%f\n" - " AliESDEvent::GetBeamType() ->%s\n" - " AliESDEvent::GetCurrentL3() ->%f\n" - " AliESDEvent::GetMagneticField()->%f\n" - " AliESDEvent::GetRunNumber() ->%d\n", - esd->GetBeamEnergy(), - esd->GetBeamType(), - esd->GetCurrentL3(), - esd->GetMagneticField(), - esd->GetRunNumber())); - - fInspector.Init(fVtxAxis); - - Print(); - fFirstEvent = false; - } - - // Some variables - UInt_t triggers; // Trigger bits - Bool_t lowFlux; // Low flux flag - UShort_t iVz; // Vertex bin from ESD - Double_t vZ; // Z coordinate from ESD - Double_t cent; // Centrality - UShort_t iVzMc; // Vertex bin from MC - Double_t vZMc; // Z coordinate of IP vertex from MC - Double_t b; // Impact parameter - Double_t cMC; // Centrality estimate from b - Int_t nPart; // Number of participants - Int_t nBin; // Number of binary collisions - Double_t phiR; // Reaction plane from MC - UShort_t nClusters;// Number of SPD clusters - // Process the data - UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, - cent, nClusters); - fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, - b, cMC, nPart, nBin, phiR); - - Bool_t isInel = triggers & AliAODForwardMult::kInel; - Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk; - - // Fill the event count histograms - if (isInel) fHEventsTr->Fill(vZMc); - if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc); - fHEvents->Fill(vZMc); - - // Now find our vertex bin object - VtxBin* bin = 0; - if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins()) - bin = static_cast(fVtxBins->At(iVzMc)); - if (!bin) { - // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc)); - return; - } + AliCentralMCCorrectionsTask::VtxBin& vb = + static_cast(bin); // Now process our input data and store in own ESD object - fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary); - + fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary); + // Get the ESD object - const AliMultiplicity* spdmult = esd->GetMultiplicity(); - + const AliMultiplicity* spdmult = esd.GetMultiplicity(); + // Count number of tracklets per bin for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) - bin->fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j)); + vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j)); //...and then the unused clusters in layer 1 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) { - Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)); - bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j)); + Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)); + vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j)); } - + // Count events - bin->fCounts->Fill(0.5); + bin.fCounts->Fill(0.5); - PostData(1, fList); + return true; } - //____________________________________________________________________ void -AliCentralMCCorrectionsTask::Terminate(Option_t*) +AliCentralMCCorrectionsTask::CreateCorrections(TList* results) { - // - // End of job - // - // Parameters: - // option Not used - // - DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output"); - - fList = dynamic_cast(GetOutputData(1)); - if (!fList) { - AliError("No output list defined"); - return; - } - - DefineBins(fList); - - // Output list - TList* output = new TList; - output->SetOwner(); - output->SetName(Form("%sResults", GetName())); - - // --- Fill correction object -------------------------------------- - AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap; - corr->SetVertexAxis(fVtxAxis); - // corr->SetEtaAxis(fEtaAxis); - - AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance; - acorr->SetVertexAxis(fVtxAxis); + fSecCorr = new AliCentralCorrSecondaryMap; + fSecCorr->SetVertexAxis(fVtxAxis); - TIter next(fVtxBins); - VtxBin* bin = 0; - UShort_t iVz = 1; - while ((bin = static_cast(next()))) - bin->Finish(fList, output, iVz++, fEffectiveCorr, corr,acorr); + fAccCorr = new AliCentralCorrAcceptance; + fAccCorr->SetVertexAxis(fVtxAxis); - output->Add(corr); - output->Add(acorr); + results->Add(fSecCorr); + results->Add(fAccCorr); +} - PostData(2, output); +//____________________________________________________________________ +Bool_t +AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin* + bin, UShort_t iVz) +{ + + AliCentralMCCorrectionsTask::VtxBin* vb = + static_cast(bin); + vb->Terminate(fList, fResults, iVz, fEffectiveCorr, + fEtaCut, fCorrCut, fSecCorr,fAccCorr); + return true; } + //____________________________________________________________________ void AliCentralMCCorrectionsTask::Print(Option_t* option) const { - std::cout << ClassName() << "\n" - << " Vertex bins: " << fVtxAxis.GetNbins() << '\n' - << " Vertex range: [" << fVtxAxis.GetXmin() - << "," << fVtxAxis.GetXmax() << "]\n" - << " Eta bins: " << fEtaAxis.GetNbins() << '\n' - << " Eta range: [" << fEtaAxis.GetXmin() - << "," << fEtaAxis.GetXmax() << "]\n" - << " # of phi bins: " << fNPhiBins + AliBaseMCCorrectionsTask::Print(option); + std::cout << " # of phi bins: " << fNPhiBins << "\n" + << " Effective corr.: " << fEffectiveCorr << "\n" + << " Eta cut-off: " << fEtaCut << "\n" + << " Acceptance cut: " << fCorrCut << std::endl; gROOT->IncreaseDirLevel(); - fInspector.Print(option); fTrackDensity.Print(option); gROOT->DecreaseDirLevel(); } //==================================================================== -const char* -AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low, - Double_t high) -{ - static TString buf; - buf = Form("vtx%+05.1f_%+05.1f", low, high); - buf.ReplaceAll("+", "p"); - buf.ReplaceAll("-", "m"); - buf.ReplaceAll(".", "d"); - return buf.Data(); -} - - -//____________________________________________________________________ AliCentralMCCorrectionsTask::VtxBin::VtxBin() - : fHits(0), - fClusters(0), - fPrimary(0), - fCounts(0) + : AliBaseMCCorrectionsTask::VtxBin(), + fHits(0), + fClusters(0) { } //____________________________________________________________________ @@ -523,21 +176,10 @@ AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low, Double_t high, const TAxis& axis, UShort_t nPhi) - : TNamed(BinName(low, high), - Form("%+5.1fcmSetXTitle("#eta"); - fPrimary->SetYTitle("#varphi [radians]"); - fPrimary->Sumw2(); - fPrimary->SetDirectory(0); - fHits = static_cast(fPrimary->Clone("hits")); fHits->SetTitle("Hits"); fHits->SetDirectory(0); @@ -545,103 +187,43 @@ AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low, fClusters = static_cast(fPrimary->Clone("clusters")); fClusters->SetTitle("Clusters"); fClusters->SetDirectory(0); - - fCounts = new TH1D("counts", "Counts", 1, 0, 1); - fCounts->SetXTitle("Events"); - fCounts->SetYTitle("# of Events"); - fCounts->SetDirectory(0); } -//____________________________________________________________________ -AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o) - : TNamed(o), - fHits(0), - fClusters(0), - fPrimary(0), - fCounts(0) -{ - if (o.fHits) { - fHits = static_cast(o.fHits->Clone()); - fHits->SetDirectory(0); - } - if (o.fClusters) { - fClusters = static_cast(o.fClusters->Clone()); - fClusters->SetDirectory(0); - } - if (o.fPrimary) { - fPrimary = static_cast(o.fPrimary->Clone()); - fPrimary->SetDirectory(0); - } - if (o.fCounts) { - fCounts = static_cast(o.fCounts->Clone()); - fCounts->SetDirectory(0); - } -} //____________________________________________________________________ -AliCentralMCCorrectionsTask::VtxBin& -AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o) +TList* +AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l) { - if (&o == this) return *this; - TNamed::operator=(o); - fHits = 0; - fPrimary = 0; - fClusters = 0; - fCounts = 0; - if (o.fHits) { - fHits = static_cast(o.fHits->Clone()); - fHits->SetDirectory(0); - } - if (o.fClusters) { - fClusters = static_cast(o.fClusters->Clone()); - fClusters->SetDirectory(0); - } - if (o.fPrimary) { - fPrimary = static_cast(o.fPrimary->Clone()); - fPrimary->SetDirectory(0); - } - if (o.fCounts) { - fCounts = static_cast(o.fCounts->Clone()); - fCounts->SetDirectory(0); - } - return *this; -} - -//____________________________________________________________________ -void -AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l) -{ - TList* d = new TList; - d->SetName(GetName()); - d->SetOwner(); - l->Add(d); + TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l); d->Add(fHits); d->Add(fClusters); - d->Add(fPrimary); - d->Add(fCounts); + + return d; } //____________________________________________________________________ void -AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input, - TList* output, - UShort_t iVz, - Bool_t effectiveCorr, - AliCentralCorrSecondaryMap* map, - AliCentralCorrAcceptance* acorr) +AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input, + TList* output, + UShort_t iVz, + Bool_t effectiveCorr, + Double_t etaCut, + Double_t accCut, + AliCentralCorrSecondaryMap* map, + AliCentralCorrAcceptance* acorr) { TList* out = new TList; out->SetName(GetName()); out->SetOwner(); output->Add(out); - + TList* l = static_cast(input->FindObject(GetName())); if (!l) { AliError(Form("List %s not found in %s", GetName(), input->GetName())); return; } - + // Get the sums TH2D* hits = static_cast(l->FindObject("hits")); TH2D* clus = static_cast(l->FindObject("clusters")); TH2D* prim = static_cast(l->FindObject("primary")); @@ -650,40 +232,108 @@ AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input, return; } - TH2D* h = 0; - if (effectiveCorr) h = static_cast(clus->Clone("bgCorr")); - else h = static_cast(hits->Clone("bgCorr")); - h->SetDirectory(0); - h->Divide(prim); + // Clone cluster and hit map + TH2D* secMapEff = static_cast(clus->Clone("secMapEff")); + TH2D* secMapHit = static_cast(hits->Clone("secMapHit")); + secMapEff->SetTitle("2^{nd} map from clusters"); + secMapEff->SetDirectory(0); + secMapHit->SetTitle("2^{nd} map from MC hits"); + secMapHit->SetDirectory(0); + + // Divide cluster and hit map with number of primaries + secMapEff->Divide(prim); + secMapHit->Divide(prim); + + // Create acceptance histograms + TH1D* accEff = new TH1D("accEff", + "Acceptance correction for SPD (from 2^{nd} map)" , + fPrimary->GetXaxis()->GetNbins(), + fPrimary->GetXaxis()->GetXmin(), + fPrimary->GetXaxis()->GetXmax()); + TH1D* accHit = static_cast(accEff->Clone("accHit")); + accHit->SetTitle("Acceptance correction for SPD (from clusters)"); + + // Diagnostics histogra, + TH2* dia = static_cast(clus->Clone("diagnostics")); + dia->SetTitle("Scaled cluster density"); + + // Total number of channels along phi and # of eta bins + Int_t nTotal = secMapHit->GetNbinsY(); + Int_t nEta = secMapHit->GetNbinsX(); + + for(Int_t xx = 1; xx <= nEta; xx++) { + Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx); + Bool_t ins = TMath::Abs(eta) <= etaCut; + Double_t mm = 0; + if (ins) { + // Find the maximum cluster signal in this phi range + for (Int_t yy = 1; yy <= nTotal; yy++) { + Double_t c = clus->GetBinContent(xx,yy); + mm = TMath::Max(mm, c); + } + } + // Count number of phi bins with enough clusters or high enough + // correction. + Int_t nOKEff = 0; + Int_t nOKHit = 0; + for(Int_t yy = 1; yy <=nTotal; yy++) { + if (!ins) { // Not inside Eta cut + secMapEff->SetBinContent(xx,yy,0.); + secMapEff->SetBinError(xx,yy,0.); + secMapHit->SetBinContent(xx,yy,0.); + secMapHit->SetBinError(xx,yy,0.); + dia->SetBinContent(xx,yy,0); + continue; + } - TH1D* acc = new TH1D(Form("SPDacc_vrtbin_%d",iVz), - "Acceptance correction for SPD" , - fPrimary->GetXaxis()->GetNbins(), - fPrimary->GetXaxis()->GetXmin(), - fPrimary->GetXaxis()->GetXmax()); - TH1F* accden = static_cast(acc->Clone(Form("%s_den",acc->GetName()))); + // Check if the background correction is large enough, or zero map + if(secMapEff->GetBinContent(xx,yy) > 0.9) { + // acc->Fill(h->GetXaxis()->GetBinCenter(xx)); + nOKEff++; + } + else { + secMapEff->SetBinContent(xx,yy,0.); + secMapEff->SetBinError(xx,yy,0.); + } - for(Int_t xx = 1; xx <=h->GetNbinsX(); xx++) { - for(Int_t yy = 1; yy <=h->GetNbinsY(); yy++) { - if(TMath::Abs(h->GetXaxis()->GetBinCenter(xx)) > 1.9) { - h->SetBinContent(xx,yy,0.); - h->SetBinError(xx,yy,0.); + // Check if the number of cluster is large enough, or zero map + Double_t c = clus->GetBinContent(xx,yy); + Double_t s = (mm < 1e-6) ? 0 : c / mm; + dia->SetBinContent(xx,yy,s); + if (s >= accCut) { + nOKHit++; } - if(h->GetBinContent(xx,yy) > 0.9) - acc->Fill(h->GetXaxis()->GetBinCenter(xx)); else { - h->SetBinContent(xx,yy,0.); - h->SetBinError(xx,yy,0.); + secMapHit->SetBinContent(xx,yy,0); + secMapHit->SetBinError(xx,yy,0); } - accden->Fill(h->GetXaxis()->GetBinCenter(xx)); } + + // Calculate acceptance as ratio of bins with enough clusters and + // total number of phi bins. + Double_t accXX = float(nOKHit) / nTotal; + if (accXX < 0.2) accXX = 0; + accHit->SetBinContent(xx, accXX); + + // Calculate acceptance as ratio of bins with large enough + // correction and total number of phi bins. + accXX = float(nOKEff) / nTotal; + if (accXX < 0.2) accXX = 0; + accEff->SetBinContent(xx, accXX); } - acc->Divide(accden); - map->SetCorrection(iVz, h); + TH2D* secMap = (effectiveCorr ? secMapEff : secMapHit); + TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff); + TH1D* acc = (effectiveCorr ? accEff : accHit); + TH1D* accAlt = (effectiveCorr ? accHit : accEff); + out->Add(secMap->Clone("secMap")); + out->Add(secMapAlt->Clone()); + out->Add(acc->Clone("acc")); + out->Add(accAlt->Clone()); + out->Add(dia->Clone()); + + map->SetCorrection(iVz, secMap); acorr->SetCorrection(iVz, acc); - out->Add(h); - out->Add(acc); } //