//
// Corrections used
#include "AliCentralMultiplicityTask.h"
+#include "AliCentralCorrectionManager.h"
+#include "AliCentralCorrAcceptance.h"
+#include "AliCentralCorrSecondaryMap.h"
#include "AliAODForwardMult.h"
#include "AliForwardUtil.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliMultiplicity.h"
#include <TROOT.h>
+#include <TSystem.h>
#include <TFile.h>
#include <TError.h>
#include <TSystem.h>
+#include <TObjArray.h>
#include <iostream>
#include <iomanip>
//====================================================================
AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
- : AliAnalysisTaskSE(name),
- fInspector("centralEventInspector"),
- fData(0),
- fList(0),
+ : AliBaseESDTask(name, "AliCentralMultiplicityTask",
+ &(AliCentralCorrectionManager::Instance())),
+ fInspector("event"),
fAODCentral(kFALSE),
- fManager(),
fUseSecondary(true),
fUseAcceptance(true),
- fFirstEventSeen(false),
- fIvz(0),
- fNClusterTracklet(0),
- fClusterPerTracklet(0),
- fNCluster(0),
- fNTracklet(0),
- fEtaMin(0),
- fEtaMax(0)
+ fIvz(0),
+ fNClusterTracklet(0),
+ fClusterPerTracklet(0),
+ fNCluster(0),
+ fNTracklet(0),
+ fVtxList(0),
+ fStore(false),
+ fHData(0)
{
//
// Constructor
//
- DefineOutput(1, TList::Class());
+ DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
fBranchNames =
"ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
"SPDVertex.,PrimaryVertex.";
}
//____________________________________________________________________
AliCentralMultiplicityTask::AliCentralMultiplicityTask()
- : AliAnalysisTaskSE(),
+ : AliBaseESDTask(),
fInspector(),
- fData(0),
- fList(0),
fAODCentral(),
- fManager(),
fUseSecondary(true),
fUseAcceptance(true),
- fFirstEventSeen(false),
- fIvz(0),
+ fIvz(0),
fNClusterTracklet(0),
fClusterPerTracklet(0),
fNCluster(0),
fNTracklet(0),
- fEtaMin(0),
- fEtaMax(0)
+ fVtxList(0),
+ fStore(false),
+ fHData(0)
{
//
// Constructor
//
+ DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
}
+
//____________________________________________________________________
-AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
- : AliAnalysisTaskSE(o),
- fInspector(o.fInspector),
- fData(o.fData),
- fList(o.fList),
- fAODCentral(o.fAODCentral),
- fManager(o.fManager),
- fUseSecondary(o.fUseSecondary),
- fUseAcceptance(o.fUseAcceptance),
- fFirstEventSeen(o.fFirstEventSeen),
- fIvz(0),
- fNClusterTracklet(o.fNClusterTracklet),
- fClusterPerTracklet(o.fClusterPerTracklet),
- fNCluster(o.fNCluster),
- fNTracklet(o.fNTracklet),
- fEtaMin(o.fEtaMin),
- fEtaMax(o.fEtaMax)
-{
- //
- // Copy constructor
- //
-}
-//____________________________________________________________________
-AliCentralMultiplicityTask&
-AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
-{
- //
- // Assignment operator
- //
- if (&o == this) return *this;
- fInspector = o.fInspector;
- fData = o.fData;
- fList = o.fList;
- fAODCentral = o.fAODCentral;
- fManager = o.fManager;
- fUseSecondary = o.fUseSecondary;
- fUseAcceptance = o.fUseAcceptance;
- fFirstEventSeen = o.fFirstEventSeen;
- fIvz = 0;
- fNClusterTracklet = o.fNClusterTracklet;
- fClusterPerTracklet= o.fClusterPerTracklet;
- fNCluster = o.fNCluster;
- fNTracklet = o.fNTracklet;
- fEtaMin = o.fEtaMin;
- fEtaMax = o.fEtaMax;
- return *this;
-}
-//____________________________________________________________________
-void AliCentralMultiplicityTask::UserCreateOutputObjects()
+void
+AliCentralMultiplicityTask::CreateBranches(AliAODHandler* ah)
{
//
// Create output objects
//
//
+ DGUARD(fDebug,1,"Create user output in AliCentralMultiplicityTask");
+
+ if (!ah)
+ //AliFatal("No AOD output handler set in analysis manager");
+ return;
- AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
- AliAODHandler* ah =
- dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
- if (!ah) AliFatal("No AOD output handler set in analysis manager");
-
-
TObject* obj = &fAODCentral;
ah->AddBranch("AliAODCentralMult", &obj);
-
-
-
- fList = new TList();
- fList->SetOwner();
-
- fInspector.DefineOutput(fList);
-
- PostData(1,fList);
+}
+//____________________________________________________________________
+Bool_t
+AliCentralMultiplicityTask::Book()
+{
+ fNeededCorrections = (AliCentralCorrectionManager::kSecondaryMap|
+ AliCentralCorrectionManager::kAcceptance);
+ return true;
}
//____________________________________________________________________
-AliESDEvent*
-AliCentralMultiplicityTask::GetESDEvent()
+Bool_t
+AliCentralMultiplicityTask::PreData(const TAxis& vertex, const TAxis& eta)
{
- //
- // Get the ESD event. IF this is the first event, initialise
- //
- AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
- if (!esd) {
- AliWarning("No ESD event found for input event");
- return 0;
- }
+ // FindEtaLimits();
+ // unsigned short s = 1;
+ TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
+ eta.GetNbins(), eta.GetXmin(), eta.GetXmax(),
+ vertex.GetNbins(),vertex.GetXmin(),
+ vertex.GetXmax());
+ hCoverage->SetDirectory(0);
+ hCoverage->SetXTitle("#eta");
+ hCoverage->SetYTitle("v_{z} [cm]");
+ hCoverage->SetZTitle("n_{bins}");
+
+ fAODCentral.Init(eta);
- // IF we've read the first event already, just return the event
- if (fFirstEventSeen) return esd;
+ UShort_t nVz = vertex.GetNbins();
+ fVtxList = new TObjArray(nVz, 1);
+ fVtxList->SetName("centMultVtxBins");
+ fVtxList->SetOwner();
- // Read the details of the rung
- fInspector.ReadRunDetails(esd);
-
- // If we weren't initialised before (i.e., in the setup), do so now.
- if (!GetManager().IsInit()) {
- GetManager().Init(fInspector.GetCollisionSystem(),
- fInspector.GetEnergy(),
- fInspector.GetField());
- AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
+ // Bool_t store = false;
+ for (Int_t v = 1; v <= nVz; v++) {
+ VtxBin* bin = new VtxBin(v, vertex.GetBinLowEdge(v),
+ vertex.GetBinUpEdge(v));
+ bin->SetupForData(fList, hCoverage, fStore);
+ fVtxList->AddAt(bin, v);
}
+ fList->Add(hCoverage);
- // Check for existence and get secondary map
- AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
- if (!secMap) AliFatal("No secondary map defined!");
- const TAxis& vaxis = secMap->GetVertexAxis();
-
- FindEtaLimits();
-
+ // Bins are
+ TArrayD bins;
+ // N-bins, loweset order, higest order, return array
+ AliForwardUtil::MakeLogScale(300, 0, 5, bins);
fNClusterTracklet = new TH2D("nClusterVsnTracklet",
"Total number of cluster vs number of tracklets",
- 100, 0, 100, 100, 0, 100);
+ bins.GetSize()-1, bins.GetArray(),
+ bins.GetSize()-1, bins.GetArray());
fNClusterTracklet->SetDirectory(0);
- fNClusterTracklet->SetXTitle("# of free clusters");
- fNClusterTracklet->SetYTitle("# of tracklets");
+ fNClusterTracklet->SetXTitle("N_{free cluster}");
+ fNClusterTracklet->SetYTitle("N_{tracklet}");
fNClusterTracklet->SetStats(0);
fList->Add(fNClusterTracklet);
fNTracklet->SetDirectory(0);
fNTracklet->Sumw2();
+ fList->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondary));
+ fList->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
+
+ fHData = static_cast<TH2D*>(fAODCentral.GetHistogram().Clone("d2Ndetadphi"));
+ fHData->SetStats(0);
+ fHData->SetDirectory(0);
+ fList->Add(fHData);
+
// Initialize the inspecto
- fInspector.Init(vaxis);
- fFirstEventSeen = kTRUE;
+ // fInspector.SetupForData(vertex);
- // Print some information
- Print();
+ fAODCentral.SetBit(AliAODCentralMult::kSecondary, fUseSecondary);
+ fAODCentral.SetBit(AliAODCentralMult::kAcceptance, fUseAcceptance);
- return esd;
-}
-//____________________________________________________________________
-void
-AliCentralMultiplicityTask::MarkEventForStore() const
-{
- // Make sure the AOD tree is filled
- AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
- AliAODHandler* ah =
- dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
- if (!ah)
- AliFatal("No AOD output handler set in analysis manager");
-
- ah->SetFillAOD(kTRUE);
+ return true;
}
-
//____________________________________________________________________
-void AliCentralMultiplicityTask::FindEtaLimits()
+Bool_t
+AliCentralMultiplicityTask::PreEvent()
{
- AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
-
- const TAxis& vaxis = secMap->GetVertexAxis();
-
- fEtaMin.Set(vaxis.GetNbins());
- fEtaMax.Set(vaxis.GetNbins());
-
- TList* hits = new TList;
- hits->SetOwner();
- hits->SetName("hitMaps");
- fList->Add(hits);
-
- TList* secs = new TList;
- secs->SetOwner();
- secs->SetName("secondaryMaps");
- fList->Add(secs);
- TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
- secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetNbins(),
- secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetXmin(),
- secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetXmax(),
- vaxis.GetNbins(),vaxis.GetXmin(),vaxis.GetXmax());
- hCoverage->SetDirectory(0);
- hCoverage->SetXTitle("#eta");
- hCoverage->SetYTitle("v_{z} [cm]");
- hCoverage->SetZTitle("n_{bins}");
- fList->Add(hCoverage);
-
- for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
- TH2D* corr = secMap->GetCorrection(UShort_t(v));
- TH1D* proj = corr->ProjectionX(Form("secCor%02d", v));
- proj->Scale(1. / corr->GetNbinsY());
- proj->SetTitle(Form("Projection of secondary correction "
- "for %+5.1f<v_{z}<%+5.1f",
- vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
- proj->SetYTitle("#LT 2^{nd} correction#GT");
- proj->SetDirectory(0);
- proj->SetMarkerStyle(20);
- proj->SetMarkerColor(kBlue+1);
- secs->Add(proj);
-
- TH2D* obg = static_cast<TH2D*>(corr->Clone(Form("secCor2DFiducial%02d",v)));
- obg->SetDirectory(0);
- secs->Add(obg);
-
- TH1D* after = static_cast<TH1D*>(proj->Clone(Form("secCorFiducial%02d",v)));
- after->SetDirectory(0);
- after->SetMarkerColor(kRed+1);
- secs->Add(after);
-
- TH2D* data = static_cast<TH2D*>(corr->Clone(Form("hitMap%02d",v)));
- //d->SetTitle(Form("hitMap%02d",v));
- data->SetTitle(Form("d^{2}N/d#eta d#phi "
- "for %+5.1f<v_{z}<%+5.1f",
- vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
- data->GetZaxis()->SetTitle("");
- data->SetMarkerColor(kBlack);
- data->SetMarkerStyle(1);
- hits->Add(data);
-
- TH1D* hAcceptance = fManager.GetAcceptanceCorrection(v);
- TH1D* accClone = static_cast<TH1D*>(hAcceptance->Clone(Form("acceptance%02d",v)));
- secs->Add(accClone);
-
- // Double_t prev = 0;
- for (Int_t e = 1; e <= proj->GetNbinsX(); e++) {
- Double_t c = proj->GetBinContent(e);
- if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
- fEtaMin[v-1] = e;
- break;
- }
- // prev = c;
- after->SetBinContent(e, 0);
- after->SetBinError(e, 0);
- for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
- obg->SetBinContent(e,nn,0);
-
- }
- for (Int_t e = proj->GetNbinsX(); e >= 1; e--) {
- Double_t c = proj->GetBinContent(e);
- if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
- fEtaMax[v-1] = e;
- break;
- }
- // prev = c;
- after->SetBinContent(e, 0);
- after->SetBinError(e, 0);
- for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
- obg->SetBinContent(e,nn,0);
-
- }
-
- for (Int_t nn = fEtaMin[v-1]; nn<=fEtaMax[v-1]; nn++) {
- hCoverage->SetBinContent(nn,v,1);
- }
-
- }
+ fAODCentral.Clear("");
+ return true;
}
-
//____________________________________________________________________
-void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
+Bool_t
+AliCentralMultiplicityTask::Event(AliESDEvent& esd)
{
//
// Process each event
// Parameters:
// option Not used
//
- fAODCentral.Clear("");
- fIvz = 0;
-
- AliESDEvent* esd = GetESDEvent();
-
+ DGUARD(fDebug,1,"Process event in AliCentralMultiplicityTask");
+ fIvz = 0;
Bool_t lowFlux = kFALSE;
UInt_t triggers = 0;
UShort_t ivz = 0;
- Double_t vz = 0;
+ TVector3 ip;
Double_t cent = -1;
UShort_t nClusters = 0;
- UInt_t found = fInspector.Process(esd, triggers, lowFlux,
- ivz, vz, cent, nClusters);
+ UInt_t found = fInspector.Process(&esd, triggers, lowFlux,
+ ivz, ip, cent, nClusters);
// No event or no trigger
- if (found & AliFMDEventInspector::kNoEvent) return;
- if (found & AliFMDEventInspector::kNoTriggers) return;
+ if (found & AliFMDEventInspector::kNoEvent) return false;
+ if (found & AliFMDEventInspector::kNoTriggers) return false;
// Make sure AOD is filled
MarkEventForStore();
- if (found == AliFMDEventInspector::kNoSPD) return;
- if (found == AliFMDEventInspector::kNoVertex) return;
- if (triggers & AliAODForwardMult::kPileUp) return;
- if (found == AliFMDEventInspector::kBadVertex) return; // Out of range
+ if (found & AliFMDEventInspector::kNoSPD) return false;
+ if (found & AliFMDEventInspector::kNoVertex) return false;
+ if (triggers & AliAODForwardMult::kPileUp) return false;
+ if (found & AliFMDEventInspector::kBadVertex) return false;
//Doing analysis
- fIvz = ivz;
- const AliMultiplicity* spdmult = esd->GetMultiplicity();
+ const AliMultiplicity* spdmult = esd.GetMultiplicity();
TH2D& aodHist = fAODCentral.GetHistogram();
+ VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivz));
+ if (!bin) return false;
+
ProcessESD(aodHist, spdmult);
- CorrectData(aodHist, ivz);
- //Producing hit maps
- TList* hitList = static_cast<TList*>(fList->FindObject("hitMaps"));
- TH2D* data = 0;
- if(hitList)
- data = static_cast<TH2D*>(hitList->At(ivz-1));
- if(data)
- data->Add(&aodHist);
+ bin->Correct(aodHist, fUseSecondary, fUseAcceptance);
- PostData(1,fList);
+ if (triggers & AliAODForwardMult::kInel)
+ fHData->Add(&(fAODCentral.GetHistogram()));
+
+ return true;
}
//____________________________________________________________________
void
AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
const AliMultiplicity* spdmult) const
{
+ DGUARD(fDebug,1,"Process the ESD in AliCentralMultiplicityTask");
fNTracklet->Reset();
fNCluster->Reset();
}
//____________________________________________________________________
-void
-AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const
-{
- // Corrections
- TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
- TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
-
- if (!hSecMap) AliFatal("No secondary map!");
- if (!hAcceptance) AliFatal("No acceptance!");
-
- if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
-
- for(Int_t nx = 1; nx <= aodHist.GetNbinsX(); nx++) {
- Float_t accCor = hAcceptance->GetBinContent(nx);
- Float_t accErr = hAcceptance->GetBinError(nx);
-
- Bool_t fiducial = true;
- if (nx < fEtaMin[vtxbin-1] || nx > fEtaMax[vtxbin-1])
- fiducial = false;
- // Bool_t etabinSeen = kFALSE;
- for(Int_t ny = 1; ny <= aodHist.GetNbinsY(); ny++) {
-#if 1
- if (!fiducial) {
- aodHist.SetBinContent(nx, ny, 0);
- aodHist.SetBinError(nx, ny, 0);
- continue;
- }
-#endif
- // Get currrent value
- Float_t aodValue = aodHist.GetBinContent(nx,ny);
- Float_t aodErr = aodHist.GetBinError(nx,ny);
-
-#if 0 // This is done once in the FindEtaBins function
- // Set underflow bin
- Float_t secCor = 0;
- if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
- if (secCor > 0.5) etabinSeen = kTRUE;
-#endif
- if (aodValue < 0.000001) {
- aodHist.SetBinContent(nx,ny, 0);
- aodHist.SetBinError(nx,ny, 0);
- continue;
- }
- if (!fUseAcceptance) continue;
-
- // Acceptance correction
- if (accCor < 0.000001) accCor = 1;
- Float_t aodNew = aodValue / accCor ;
- Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
- TMath::Power(accErr/accCor,2) );
- aodHist.SetBinContent(nx,ny, aodNew);
- //test
- aodHist.SetBinError(nx,ny,error);
- aodHist.SetBinError(nx,ny,aodErr);
- }
- //Filling underflow bin if we eta bin is in range
- if (fiducial) aodHist.SetBinContent(nx,0, 1.);
- // if (etabinSeen) aodHist.SetBinContent(nx,0, 1.);
- }
-}
-
-//____________________________________________________________________
-void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
+Bool_t
+AliCentralMultiplicityTask::Finalize()
{
//
// End of job
// Parameters:
// option Not used
//
+ DGUARD(fDebug,1,"Process merged output in AliCentralMultiplicityTask");
+
+ Double_t nTr = 0, nTrVtx = 0, nAcc = 0;
+ MakeSimpledNdeta(fList, fResults, nTr, nTrVtx, nAcc);
+ AliInfoF("\n"
+ "\t# events w/trigger: %f\n"
+ "\t# events w/trigger+vertex: %f\n"
+ "\t# events accepted y cuts: %f",
+ nTr, nTrVtx, nAcc);
+ return true;
}
+
+//____________________________________________________________________
+Bool_t
+AliCentralMultiplicityTask::MakeSimpledNdeta(const TList* input,
+ TList* output,
+ Double_t& nTr,
+ Double_t& nTrVtx,
+ Double_t& nAcc)
+{
+ // Get our histograms from the container
+ TH1I* hEventsTr = 0;
+ TH1I* hEventsTrVtx = 0;
+ TH1I* hEventsAcc = 0;
+ TH1I* hTriggers = 0;
+ if (!GetEventInspector().FetchHistograms(input,
+ hEventsTr,
+ hEventsTrVtx,
+ hEventsAcc,
+ hTriggers)) {
+ AliError(Form("Didn't get histograms from event selector "
+ "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p,hTriggers=%p)",
+ hEventsTr, hEventsTrVtx, hEventsAcc, hTriggers));
+ input->ls();
+ return false;
+ }
+ nTr = hEventsTr->Integral();
+ nTrVtx = hEventsTrVtx->Integral();
+ nAcc = hEventsAcc->Integral();
+ Double_t vtxEff = nTrVtx / nTr;
+ TH2D* hData = static_cast<TH2D*>(input->FindObject("d2Ndetadphi"));
+ if (!hData) {
+ AliError(Form("Couldn't get our summed histogram from output "
+ "list %s (d2Ndetadphi=%p)", input->GetName(), hData));
+ input->ls();
+ return false;
+ }
+
+ Int_t nY = hData->GetNbinsY();
+ TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, nY, "e");
+ TH1D* dNdeta_ = hData->ProjectionX("dNdeta_", 1, nY, "e");
+ TH1D* norm = hData->ProjectionX("norm", 0, 0, "");
+ TH1D* phi = hData->ProjectionX("phi", nY+1, nY+1, "");
+ dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
+ dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ dNdeta->SetMarkerColor(kRed+1);
+ dNdeta->SetMarkerStyle(20);
+ dNdeta->SetDirectory(0);
+
+ dNdeta_->SetTitle("dN_{ch}/d#eta in the forward regions");
+ dNdeta_->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
+ dNdeta_->SetMarkerColor(kMagenta+1);
+ dNdeta_->SetMarkerStyle(21);
+ dNdeta_->SetDirectory(0);
+
+ norm->SetTitle("Normalization to #eta coverage");
+ norm->SetYTitle("#eta coverage");
+ norm->SetLineColor(kBlue+1);
+ norm->SetMarkerColor(kBlue+1);
+ norm->SetMarkerStyle(21);
+ norm->SetFillColor(kBlue+1);
+ norm->SetFillStyle(3005);
+ norm->SetDirectory(0);
+
+ phi->SetTitle("Normalization to #phi acceptance");
+ phi->SetYTitle("#phi acceptance");
+ phi->SetLineColor(kGreen+1);
+ phi->SetMarkerColor(kGreen+1);
+ phi->SetMarkerStyle(20);
+ phi->SetFillColor(kGreen+1);
+ phi->SetFillStyle(3004);
+ // phi->Scale(1. / nAcc);
+ phi->SetDirectory(0);
+
+ // dNdeta->Divide(norm);
+ dNdeta->Divide(phi);
+ dNdeta->SetStats(0);
+ dNdeta->Scale(vtxEff, "width");
+
+ dNdeta_->Divide(norm);
+ dNdeta_->SetStats(0);
+ dNdeta_->Scale(vtxEff, "width");
+
+ output->Add(dNdeta);
+ output->Add(dNdeta_);
+ output->Add(norm);
+ output->Add(phi);
+
+ return true;
+}
+
+#define PF(N,V,...) \
+ AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
+#define PFB(N,FLAG) \
+ do { \
+ AliForwardUtil::PrintName(N); \
+ std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
+ } while(false)
+#define PFV(N,VALUE) \
+ do { \
+ AliForwardUtil::PrintName(N); \
+ std::cout << (VALUE) << std::endl; } while(false)
+
//____________________________________________________________________
void
AliCentralMultiplicityTask::Print(Option_t* option) const
// Parameters:
// option Not used
//
- std::cout << ClassName() << ": " << GetName() << "\n"
- << std::boolalpha
- << " Use secondary correction: " << fUseSecondary << '\n'
- << " Use acceptance correction: " << fUseAcceptance << '\n'
- << " Off-line trigger mask: 0x"
- << std::hex << std::setfill('0')
- << std::setw (8) << fOfflineTriggerMask
- << std::dec << std::setfill (' ')
- << std::noboolalpha << std::endl;
- AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
- if (secMap) {
- const TAxis& vaxis = secMap->GetVertexAxis();
- std::cout << " Eta ranges:\n"
- << " Vertex | Eta bins\n"
- << " bin range | \n"
- << " ----------------+-----------" << std::endl;
- for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
- std::cout << " " << std::setw(2) << v << " "
- << std::setw(5) << vaxis.GetBinLowEdge(v) << "-"
- << std::setw(5) << vaxis.GetBinUpEdge(v) << " | "
- << std::setw(3) << fEtaMin[v-1] << "-"
- << std::setw(3) << fEtaMax[v-1] << std::endl;
+ AliBaseESDTask::Print(option);
+ PFB("Use secondary correction", fUseSecondary);
+ PFB("Use acceptance correction", fUseAcceptance);
+
+ AliCentralCorrectionManager& ccm =
+ AliCentralCorrectionManager::Instance();
+ if (ccm.IsInit()) {
+ const AliCentralCorrSecondaryMap* secMap = ccm.GetSecondaryMap();
+ if (secMap) {
+ const TAxis& vaxis = secMap->GetVertexAxis();
+ // fVtxList->ls();
+ std::cout << " Eta ranges:\n"
+ << " Vertex | Eta bins\n"
+ << " bin range | \n"
+ << " ----------------+-----------" << std::endl;
+ for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
+ VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(v));
+ if (!bin) continue;
+ bin->Print();
+ }
}
}
gROOT->IncreaseDirLevel();
- fManager.Print(option);
- fInspector.Print(option);
+ ccm.Print(option);
+ // fInspector.Print(option);
gROOT->DecreaseDirLevel();
}
+
//====================================================================
-AliCentralMultiplicityTask::Manager::Manager() :
- fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
- fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
- fAcceptance(),
- fSecmap(),
- fAcceptanceName("centralacceptance"),
- fSecMapName("centralsecmap"),
- fIsInit(kFALSE)
+AliCentralMultiplicityTask::VtxBin::VtxBin(Int_t iVz,
+ Double_t minIpZ,
+ Double_t maxIpZ)
+ : fId(iVz),
+ fMinIpZ(minIpZ),
+ fMaxIpZ(maxIpZ),
+ fEtaMin(999),
+ fEtaMax(0),
+ fSec(0),
+ fAcc(0),
+ fHits(0)
{
- //
- // Constructor
- //
}
//____________________________________________________________________
-AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
- :fAcceptancePath(o.fAcceptancePath),
- fSecMapPath(o.fSecMapPath),
- fAcceptance(o.fAcceptance),
- fSecmap(o.fSecmap),
- fAcceptanceName(o.fAcceptanceName),
- fSecMapName(o.fSecMapName),
- fIsInit(o.fIsInit)
+AliCentralMultiplicityTask::VtxBin::VtxBin(const VtxBin& o)
+ : TObject(o),
+ fId(o.fId),
+ fMinIpZ(o.fMinIpZ),
+ fMaxIpZ(o.fMaxIpZ),
+ fEtaMin(o.fEtaMin),
+ fEtaMax(o.fEtaMax),
+ fSec(o.fSec),
+ fAcc(o.fAcc),
+ fHits(o.fHits)
{
- //
- // Copy Constructor
- //
}
//____________________________________________________________________
-AliCentralMultiplicityTask::Manager&
-AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
+AliCentralMultiplicityTask::VtxBin&
+AliCentralMultiplicityTask::VtxBin::operator=(const VtxBin& o)
{
- //
- // Assignment operator
- //
- if (&o == this) return *this;
- fAcceptancePath = o.fAcceptancePath;
- fSecMapPath = o.fSecMapPath;
- fAcceptance = o.fAcceptance;
- fSecmap = o.fSecmap;
- fAcceptanceName = o.fAcceptanceName;
- fSecMapName = o.fSecMapName;
- fIsInit = o.fIsInit;
+ if (&o == this) return *this;
+ fId = o.fId;
+ fMinIpZ = o.fMinIpZ;
+ fMaxIpZ = o.fMaxIpZ;
+ fEtaMin = o.fEtaMin;
+ fEtaMax = o.fEtaMax;
+ fSec = o.fSec;
+ fAcc = o.fAcc;
+ fHits = o.fHits;
+
return *this;
}
//____________________________________________________________________
-const char*
-AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
- UShort_t sys,
- UShort_t sNN,
- Short_t field) const
+const char*
+AliCentralMultiplicityTask::VtxBin::GetName() const
{
- //
- // Get full path name to object file
- //
- // Parameters:
- // what What to get
- // sys Collision system
- // sNN Center of mass energy
- // field Magnetic field
- //
- // Return:
- //
- //
- return Form("%s/%s",
- what == 0 ? GetSecMapPath() : GetAcceptancePath(),
- GetFileName(what, sys, sNN, field));
+ return Form("%c%03d_%c%03d",
+ (fMinIpZ >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fMinIpZ)),
+ (fMaxIpZ >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fMaxIpZ)));
}
//____________________________________________________________________
-const char*
-AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
- UShort_t sys,
- UShort_t sNN,
- Short_t field) const
+void
+AliCentralMultiplicityTask::VtxBin::SetupForData(TList* l,
+ TH2* coverage,
+ Bool_t store)
{
- //
- // Get the full path name
- //
- // Parameters:
- // what What to get
- // sys Collision system
- // sNN Center of mass energy
- // field Magnetic field
- //
- // Return:
- //
- //
- // Must be static - otherwise the data may disappear on return from
- // this member function
- static TString fname = "";
-
- switch(what) {
- case 0: fname = fSecMapName; break;
- case 1: fname = fAcceptanceName; break;
- default:
- ::Error("GetFileName",
- "Invalid indentifier %d for central object, must be 0 or 1!", what);
- break;
+ TList* out = 0;
+ if (store) {
+ out = new TList;
+ out->SetName(GetName());
+ out->SetOwner();
+ l->Add(out);
}
- fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
- AliForwardUtil::CollisionSystemString(sys),
- sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
-
- return fname.Data();
-}
-//____________________________________________________________________
-TH2D*
-AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
-{
- //
- // Get the secondary map
- //
- // Parameters:
- // vtxbin
- //
- // Return:
- //
- //
- if (!fSecmap) {
- ::Warning("GetSecMapCorrection","No secondary map defined");
- return 0;
- }
- return fSecmap->GetCorrection(vtxbin);
-}
-//____________________________________________________________________
-TH1D*
-AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
- const
-{
- //
- // Get the acceptance correction
- //
- // Parameters:
- // vtxbin
- //
- // Return:
- //
- //
- if (!fAcceptance) {
- ::Warning("GetAcceptanceCorrection","No acceptance map defined");
- return 0;
- }
- return fAcceptance->GetCorrection(vtxbin);
-}
+ AliCentralCorrectionManager& ccm =
+ AliCentralCorrectionManager::Instance();
-//____________________________________________________________________
-void
-AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
- UShort_t sNN,
- Short_t field)
-{
- //
- // Initialize
- //
- // Parameters:
- // sys Collision system (1: pp, 2: PbPb, 3: pPb)
- // sNN Center of mass energy per nucleon pair [GeV]
- // field Magnetic field [kG]
- //
- if(fIsInit) ::Warning("Init","Already initialised - overriding...");
-
- TFile fsec(GetFullFileName(0,sys,sNN,field));
- fSecmap =
- dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
- if(!fSecmap) {
- ::Error("Init", "no central Secondary Map found!") ;
- return;
+ // Clean-up
+ if (fSec) {
+ // delete fSec;
+ fSec = 0;
}
- TFile facc(GetFullFileName(1,sys,sNN,field));
- fAcceptance =
- dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
- if(!fAcceptance) {
- ::Error("Init", "no central Acceptance found!") ;
- return;
+ if (fAcc) {
+ // delete fAcc;
+ fAcc = 0;
}
-
- if(fSecmap && fAcceptance) {
- fIsInit = kTRUE;
- ::Info("Init",
- "Central Manager initialised for %s, energy %dGeV, field %dkG",
- sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "unknown", sNN,field);
- }
-}
-//____________________________________________________________________
-Bool_t
-AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what,
- UShort_t sys,
- UShort_t sNN,
- Short_t fld,
- TObject* obj,
- Bool_t full) const
-{
- //
- // Write correction output to (a temporary) file
- //
- // Parameters:
- // What What to write
- // sys Collision system (1: pp, 2: PbPb, 3: pPb)
- // sNN Center of mass energy per nucleon (GeV)
- // fld Field (kG)
- // obj Object to write
- // full if true, write to full path, otherwise locally
- //
- // Return:
- // true on success.
- TString ofName;
- if (!full)
- ofName = GetFileName(what, sys, sNN, fld);
- else
- ofName = GetFullFileName(what, sys, sNN, fld);
- if (ofName.IsNull()) {
- AliErrorGeneral("Manager",Form("Unknown object type %d", what));
- return false;
+ // Get secondary correction and make a projection onto eta
+ TH2* sec = ccm.GetSecondaryMap()->GetCorrection(UShort_t(fId));
+ TH1* acc = ccm.GetAcceptance()->GetCorrection(UShort_t(fId));
+ fSec = static_cast<TH2*>(sec->Clone());
+ fAcc = static_cast<TH1*>(acc->Clone());
+ fSec->SetDirectory(0);
+ fAcc->SetDirectory(0);
+
+ TH1D* proj = fSec->ProjectionX("secondary");
+ proj->SetDirectory(0);
+ proj->Scale(1. / fSec->GetNbinsY());
+
+ // Find lower bound on eta
+ fEtaMin = proj->GetNbinsX();
+ for (Int_t e = 1; e <= proj->GetNbinsX(); e++) {
+ Double_t c = proj->GetBinContent(e);
+ if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
+ fEtaMin = e;
+ break;
+ }
}
- TFile* output = TFile::Open(ofName, "RECREATE");
- if (!output) {
- AliErrorGeneral("Manager",Form("Failed to open file %s", ofName.Data()));
- return false;
+ // Find upper bound on eta
+ fEtaMax = 1;
+ for (Int_t e = proj->GetNbinsX(); e >= 1; e--) {
+ Double_t c = proj->GetBinContent(e);
+ if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
+ fEtaMax = e;
+ break;
+ }
+ }
+ // Fill our coverage histogram
+ for (Int_t nn = fEtaMin; nn<=fEtaMax; nn++) {
+ coverage->SetBinContent(nn,fId,1);
}
- TString oName(GetObjectName(what));
- Int_t ret = obj->Write(oName);
- if (ret <= 0) {
- AliErrorGeneral("Manager",Form("Failed to write %p to %s/%s (%d)",
- obj, ofName.Data(), oName.Data(), ret));
- return false;
+ if (!store) {
+ // If we're not asked to store anything, clean-up, and get out
+ delete proj;
+ return;
}
- ret = output->Write();
- if (ret < 0) {
- AliErrorGeneral("Manager",
- Form("Failed to write %s to disk (%d)", ofName.Data(),ret));
- return false;
+ // Modify the title of the projection
+ proj->SetTitle(Form("Projection of secondary correction "
+ "for %+5.1f<v_{z}<%+5.1f",fMinIpZ, fMaxIpZ));
+ proj->SetYTitle("#LT 2^{nd} correction#GT");
+ proj->SetMarkerStyle(20);
+ proj->SetMarkerColor(kBlue+1);
+ out->Add(proj);
+
+ // Make some histograms to store diagnostics
+ TH2D* obg = static_cast<TH2D*>(fSec->Clone("secondaryMapFiducial"));
+ obg->SetTitle(Form("%s - fiducial volume", obg->GetTitle()));
+ obg->GetYaxis()->SetTitle("#varphi");
+ obg->SetDirectory(0);
+ out->Add(obg);
+
+ TH1D* after = static_cast<TH1D*>(proj->Clone("secondaryFiducial"));
+ after->SetDirectory(0);
+ after->GetYaxis()->SetTitle("#LT 2^{nd} correction#GT");
+ after->SetTitle(Form("%s - fiducial volume", after->GetTitle()));
+ after->SetMarkerColor(kRed+1);
+ out->Add(after);
+
+ if (fHits) {
+ // delete fHits;
+ fHits = 0;
+ }
+ fHits = static_cast<TH2D*>(fSec->Clone("hitMap"));
+ fHits->SetDirectory(0);
+ fHits->SetTitle(Form("d^{2}N/d#eta d#phi for %+5.1f<v_{z}<%+5.1f",
+ fMinIpZ, fMaxIpZ));
+ fHits->GetYaxis()->SetTitle("#varphi");
+ fHits->GetZaxis()->SetTitle("d^{2}N/d#eta d#varphi");
+ fHits->SetMarkerColor(kBlack);
+ fHits->SetMarkerStyle(1);
+ out->Add(fHits);
+
+ // Get the acceptance, and store that
+ TH1D* accClone = static_cast<TH1D*>(fAcc->Clone("acceptance"));
+ accClone->SetTitle(Form("Acceptance for %+5.1f<v_{z}<%+5.1f",
+ fMinIpZ, fMaxIpZ));
+ accClone->SetDirectory(0);
+ out->Add(accClone);
+
+ // Now zero content outside our eta range
+ for (Int_t e = 1; e < fEtaMin; e++) {
+ after->SetBinContent(e, 0);
+ after->SetBinError(e, 0);
+ for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
+ obg->SetBinContent(e,nn,0);
}
- output->ls();
- output->Close();
-
- TString cName(obj->IsA()->GetName());
- AliInfoGeneral("Manager",
- Form("Wrote %s object %s to %s\n",
- cName.Data(),oName.Data(), ofName.Data()));
- if (!full) {
- TString dName(GetFileDir(what));
- AliInfoGeneral("Manager",
- Form("%s should be copied to %s\n"
- "Do for example\n\t"
- "aliroot $ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/"
- "MoveCorrections.C\\(%d\\)\nor\n\t"
- "cp %s %s/",
- ofName.Data(),dName.Data(),
- what, ofName.Data(),
- gSystem->ExpandPathName(dName.Data())));
-
+ for (Int_t e = fEtaMax+1; e <= proj->GetNbinsX(); e++) {
+ after->SetBinContent(e, 0);
+ after->SetBinError(e, 0);
+ for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
+ obg->SetBinContent(e,nn,0);
}
- return true;
}
+//____________________________________________________________________
+void
+AliCentralMultiplicityTask::VtxBin::Correct(TH2D& aodHist,
+ Bool_t useSecondary,
+ Bool_t useAcceptance,
+ Bool_t sum) const
+{
+ if (useSecondary && fSec) aodHist.Divide(fSec);
+
+ Int_t nY = aodHist.GetNbinsY();
+ for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) {
+ Bool_t fiducial = true;
+ if (ix < fEtaMin || ix > fEtaMax) fiducial = false;
+ // Bool_t etabinSeen = kFALSE;
+
+ Double_t eta = aodHist.GetXaxis()->GetBinCenter(ix);
+ Int_t iax = fAcc->GetXaxis()->FindBin(eta);
+ Float_t accCor = fAcc->GetBinContent(iax);
+ // For test
+ // Float_t accErr = fAcc->GetBinError(ix);
+
+ // Loop over phi
+ for(Int_t iy = 1; iy <= nY; iy++) {
+ // If outside our fiducial volume, zero content
+ if (!fiducial) {
+ aodHist.SetBinContent(ix, iy, 0);
+ aodHist.SetBinError(ix, iy, 0);
+ continue;
+ }
+ // Get currrent value
+ Float_t aodValue = aodHist.GetBinContent(ix,iy);
+ Float_t aodErr = aodHist.GetBinError(ix,iy);
+ // Ignore very small values
+ if (aodValue < 0.000001) {
+ aodHist.SetBinContent(ix,iy, 0);
+ aodHist.SetBinError(ix,iy, 0);
+ continue;
+ }
+ if (!useAcceptance) continue;
+
+ // Acceptance correction
+ Float_t accTmp = accCor;
+ if (accTmp < 0.000001) accTmp = 1;
+ Float_t aodNew = aodValue / accTmp ;
+ aodHist.SetBinContent(ix,iy, aodNew);
+ aodHist.SetBinError(ix,iy,aodErr);
+ // - Test -
+ // Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
+ // TMath::Power(accErr/accCor,2) );
+ // test - aodHist.SetBinError(ix,iy,error);
+ } // for (iy)
+ //Filling underflow bin if we eta bin is in range
+ if (fiducial) {
+ aodHist.SetBinContent(ix,0, 1.);
+ aodHist.SetBinContent(ix,nY+1, accCor);
+ }
+ } // for (ix)
+ if (sum && fHits) fHits->Add(&aodHist);
+}
+
//____________________________________________________________________
-void
-AliCentralMultiplicityTask::Manager::Print(Option_t* option) const
+void
+AliCentralMultiplicityTask::VtxBin::Print(Option_t* /*option*/) const
{
- //
- // Print information to standard output
- //
- std::cout << " AliCentralMultiplicityTask::Manager\n"
- << std::boolalpha
- << " Initialized: " << fIsInit << '\n'
- << " Acceptance path: " << fAcceptancePath << '\n'
- << " Acceptance name: " << fAcceptanceName << '\n'
- << " Acceptance: " << fAcceptance << '\n'
- << " Secondary path: " << fSecMapPath << '\n'
- << " Secondary name: " << fSecMapName << '\n'
- << " Secondary map: " << fSecmap
- << std::noboolalpha << std::endl;
- if (fAcceptance) fAcceptance->Print(option);
- if (fSecmap) fSecmap->Print(option);
+ std::cout << " "
+ << std::setw(2) << fId << " "
+ << std::setw(5) << fMinIpZ << "-"
+ << std::setw(5) << fMaxIpZ << " | "
+ << std::setw(3) << fEtaMin << "-"
+ << std::setw(3) << fEtaMax << std::endl;
}
//
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff