[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwarddNdetaTask.cxx

//====================================================================
#include "AliForwarddNdetaTask.h"
#include <TMath.h>
#include <TH2D.h>
#include <TH1D.h>
#include <THStack.h>
#include <TList.h>
#include <AliAnalysisManager.h>
#include <AliAODEvent.h>
#include <AliAODHandler.h>
#include <AliAODInputHandler.h>
#include "AliForwardUtil.h"
#include <AliAODForwardMult.h>

//____________________________________________________________________
AliForwarddNdetaTask::AliForwarddNdetaTask()
  : AliAnalysisTaskSE(), 
    fSumForward(0),	//  Sum of histograms 
    fSumForwardMC(0),	//  Sum of MC histograms (if any)
    fSumPrimary(0),	//  Sum of primary histograms
    fSumCentral(0),	//  Sum of central histograms
    fCentral(0),	// Cache of central histogram
    fPrimary(0),	// Cache of primary histogram
    fSums(0),		// Container of sums 
    fOutput(0),		// Container of outputs 
    fTriggers(0),	// Histogram of triggers 
    fVtxMin(0),		// Minimum v_z
    fVtxMax(0),		// Maximum v_z
    fTriggerMask(0),    // Trigger mask 
    fRebin(0),       	// Rebinning factor 
    fCutEdges(false),
    fSNNString(0),
    fSysString(0)
{}

//____________________________________________________________________
AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
  : AliAnalysisTaskSE("Forward"), 
    fSumForward(0),	//  Sum of histograms 
    fSumForwardMC(0),	//  Sum of MC histograms (if any)
    fSumPrimary(0),	//  Sum of primary histograms
    fSumCentral(0),	//  Sum of central histograms
    fCentral(0),	// Cache of central histogram
    fPrimary(0),	// Cache of primary histogram
    fSums(0),		// Container of sums 
    fOutput(0),		// Container of outputs 
    fTriggers(0),	// Histogram of triggers 
    fVtxMin(-10),	// Minimum v_z
    fVtxMax(10),	// Maximum v_z
    fTriggerMask(AliAODForwardMult::kInel), 
    fRebin(5),		// Rebinning factor 
    fCutEdges(false),
    fSNNString(0),
    fSysString(0)
{
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class()); 
  DefineOutput(2, TList::Class()); 
}

//____________________________________________________________________
AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
  : AliAnalysisTaskSE(o), 
    fSumForward(o.fSumForward),	   // Sum of histograms 
    fSumForwardMC(o.fSumForwardMC),// Sum of MC histograms (if any)
    fSumPrimary(o.fSumPrimary),	   // Sum of primary histograms
    fSumCentral(o.fSumCentral),	   // Sum of central histograms
    fCentral(o.fCentral),          // Cache of central histogram
    fPrimary(o.fPrimary),          // Cache of primary histogram
    fSums(o.fSums),		   // Container of sums 
    fOutput(o.fOutput),		   // Container of outputs 
    fTriggers(o.fTriggers),	   // Histogram of triggers 
    fVtxMin(o.fVtxMin),		   // Minimum v_z
    fVtxMax(o.fVtxMax),		   // Maximum v_z
    fTriggerMask(o.fTriggerMask),  // Trigger mask 
    fRebin(o.fRebin),		   // Rebinning factor 
    fCutEdges(o.fCutEdges),	   // Whether to cut edges when rebinning
    fSNNString(o.fSNNString),	   //  
    fSysString(o.fSysString)	   //  
{}

//____________________________________________________________________
AliForwarddNdetaTask::~AliForwarddNdetaTask()
{
  if (fSums) { 
    fSums->Delete();
    delete fSums;
    fSums = 0;
  }
  if (fOutputs) { 
    fOutputs->Delete();
    delete fOutputs;
    fOutputs = 0;
  }
  if (fCentral)    delete fCentral;
  if (fPrimary)    delete fPrimary;
}

//________________________________________________________________________
void 
AliForwarddNdetaTask::SetTriggerMask(const char* mask)
{
  UShort_t    trgMask = 0;
  TString     trgs(mask);
  trgs.ToUpper();
  TObjString* trg;
  TIter       next(trgs.Tokenize(" ,|"));
  while ((trg = static_cast<TObjString*>(next()))) { 
    TString s(trg->GetString());
    if      (s.IsNull()) continue;
    if      (s.CompareTo("INEL")  == 0) trgMask = AliAODForwardMult::kInel;
    else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
    else if (s.CompareTo("NSD")   == 0) trgMask = AliAODForwardMult::kNSD;
    else 
      Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
  }
  if (trgMask == 0) trgMask = 1;
  SetTriggerMask(trgMask);
}

//________________________________________________________________________
void 
AliForwarddNdetaTask::UserCreateOutputObjects()
{
  // Create histograms
  // Called once (on the worker node)

  fOutput = new TList;
  fOutput->SetName(Form("%s_result", GetName()));
  fOutput->SetOwner();

  fSums = new TList;
  fSums->SetName(Form("%s_sums", GetName()));
  fSums->SetOwner();


  fTriggers = new TH1D("triggers", "Number of triggers", 
		       kAccepted, 1, kAccepted);
  fTriggers->SetYTitle("# of events");
  fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
  fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
  fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
  fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
  fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
  fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
  fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
  fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
  fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
  fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
  fTriggers->SetFillColor(kRed+1);
  fTriggers->SetFillStyle(3001);
  fTriggers->SetStats(0);
  fSums->Add(fTriggers);

  fSNNString = new TNamed("sNN", "");
  fSysString = new TNamed("sys", "");
  fSums->Add(fSNNString);
  fSums->Add(fSysString);

  // Check that we have an AOD input handler 
  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODInputHandler* ah = 
    dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
  if (!ah) AliFatal("No AOD input handler set in analysis manager");

  // Post data for ALL output slots >0 here, to get at least an empty histogram
  PostData(1, fSums); 
}

//____________________________________________________________________
TH2D*
AliForwarddNdetaTask::CloneHist(TH2D* in, const char* name) 
{
  if (!in) return 0;
  TH2D* ret = static_cast<TH2D*>(in->Clone(name));
  ret->SetDirectory(0);
  ret->Sumw2();
  fSums->Add(ret);

  return ret;
}

//____________________________________________________________________
void 
AliForwarddNdetaTask::UserExec(Option_t *) 
{
  // Main loop
  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
  if (!aod) {
    AliError("Cannot get the AOD event");
    return;
  }  

  // Get objects from the event structure 
  TObject* oForward   = aod->FindListObject("Forward");
  TObject* oForwardMC = aod->FindListObject("ForwardMC");
  TObject* oPrimary   = aod->FindListObject("primary");
  TObject* oCentral   = aod->FindListObject("Central");

  // We should have a forward object at least 
  if (!oForward) { 
    AliWarning("No Forward object found AOD");
    return;
  }

  // Cast to good types 
  AliAODForwardMult* forward   = static_cast<AliAODForwardMult*>(oForward);
  AliAODForwardMult* forwardMC = static_cast<AliAODForwardMult*>(oForwardMC);
  TH2D*              primary   = static_cast<TH2D*>(oPrimary);
  TH2D*              central   = static_cast<TH2D*>(oCentral);

  static bool first = true;
  if (first) { 
    UShort_t sNN = forward->GetSNN();
    UShort_t sys = forward->GetSystem();
    fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
    fSNNString->SetUniqueID(sNN);
    fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
    fSysString->SetUniqueID(sys);
    first = false;
  }

  // Create our sum histograms 
  if (!fSumForward) fSumForward = CloneHist(&forward->GetHistogram(),"forward");
  if (!fSumForwardMC && forwardMC) 
    fSumForwardMC = CloneHist(&forwardMC->GetHistogram(),"forwardMC");
  if (!fSumPrimary && primary) fSumPrimary = CloneHist(primary, "primary");
  if (!fSumCentral && central) fSumCentral = CloneHist(central, "central");

  // Add contribtion from MC 
  if (primary) fSumPrimary->Add(primary);
  
  // Count event 
  fTriggers->AddBinContent(kAll);
  if (forward->IsTriggerBits(AliAODForwardMult::kB)) 
    fTriggers->AddBinContent(kB);
  if (forward->IsTriggerBits(AliAODForwardMult::kA)) 
    fTriggers->AddBinContent(kA);
  if (forward->IsTriggerBits(AliAODForwardMult::kC)) 
    fTriggers->AddBinContent(kC);
  if (forward->IsTriggerBits(AliAODForwardMult::kE)) 
    fTriggers->AddBinContent(kE);
  if (forward->IsTriggerBits(AliAODForwardMult::kInel)) 
    fTriggers->AddBinContent(kMB);

  // Check if we have an event of interest. 
  if (!forward->IsTriggerBits(fTriggerMask)) return;
  fTriggers->AddBinContent(kWithTrigger);

  // Check that we have a valid vertex
  if (!forward->HasIpZ()) return;
  fTriggers->AddBinContent(kWithVertex);

  // Check that vertex is within cuts 
  if (!forward->InRange(fVtxMin, fVtxMax)) return;
  fTriggers->AddBinContent(kAccepted);

  // Add contribution 
  fSumForward->Add(&(forward->GetHistogram()));
  if (fSumForwardMC) fSumForwardMC->Add(&(forwardMC->GetHistogram()));
  if (fSumPrimary)   fSumPrimary->Add(primary);
  if (fSumCentral)   fSumCentral->Add(central);

  PostData(1, fSums);
}

//________________________________________________________________________
void 
AliForwarddNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
					  const char* title, const char* ytitle)
{
  h->SetTitle(title);
  h->SetMarkerColor(colour);
  h->SetMarkerStyle(marker);
  h->SetMarkerSize(1);
  h->SetFillStyle(0);
  h->SetYTitle(ytitle);
  h->GetXaxis()->SetTitleFont(132);
  h->GetXaxis()->SetLabelFont(132);
  h->GetXaxis()->SetNdivisions(10);
  h->GetYaxis()->SetTitleFont(132);
  h->GetYaxis()->SetLabelFont(132);
  h->GetYaxis()->SetNdivisions(10);
  h->GetYaxis()->SetDecimals();
  h->SetStats(0);
}
  
//________________________________________________________________________
void 
AliForwarddNdetaTask::Terminate(Option_t *) 
{
  // Draw result to screen, or perform fitting, normalizations Called
  // once at the end of the query
        
  fSums = dynamic_cast<TList*> (GetOutputData(1));
  if(!fSums) {
    AliError("Could not retrieve TList fSums"); 
    return; 
  }
  
  if (!fOutput) { 
    fOutput = new TList;
    fOutput->SetName(Form("%s_result", GetName()));
    fOutput->SetOwner();
  }

  fSumForward     = static_cast<TH2D*>(fSums->FindObject("forward"));
  fSumForwardMC   = static_cast<TH2D*>(fSums->FindObject("forwardMC"));
  fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("primary"));
  fSumCentral     = static_cast<TH2D*>(fSums->FindObject("central"));
  fTriggers       = static_cast<TH1D*>(fSums->FindObject("triggers"));
  fSNNString      = static_cast<TNamed*>(fSums->FindObject("sNN"));
  fSysString      = static_cast<TNamed*>(fSums->FindObject("sys"));

  if (!fTriggers) { 
    AliError("Couldn't find histogram 'triggers' in list");
    return;
  }
  if (!fSumForward) { 
    AliError("Couldn't find histogram 'forward' in list");
    return;
  }

  Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
  Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
  Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
  Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
  Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
  Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
  Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
  Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
  Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
  Int_t nGood       = nB - nA - nC + 2 * nE;
  Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
  AliInfo(Form("Total of %9d events\n"
	       "                   of these %9d are minimum bias\n"
	       "                   of these %9d has a %s trigger\n" 
	       "                   of these %9d has a vertex\n" 
	       "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
	       "                   Triggers by type:\n"
	       "                     B   = %9d\n"
	       "                     A|C = %9d (%9d+%-9d)\n"
	       "                     E   = %9d\n"
	       "                   Implies %9d good triggers\n"
	       "                   Vertex efficiency: %f",
	       nAll, nMB, nTriggered, 
	       AliAODForwardMult::GetTriggerString(fTriggerMask),
	       nWithVertex, nAccepted,
	       fVtxMin, fVtxMax, 
	       nB, nA+nC, nA, nC, nE, nGood, vtxEff));
  
  // Get acceptance normalisation from underflow bins 
  TH1D* normForward   = fSumForward->ProjectionX("normForward", 0, 1, "");
  // Project onto eta axis - _ignoring_underflow_bins_!
  TH1D* dndetaForward = fSumForward->ProjectionX("dndetaForward", 1, -1, "e");
  // Normalize to the acceptance 
  dndetaForward->Divide(normForward);
  // Scale by the vertex efficiency 
  dndetaForward->Scale(vtxEff, "width");

  SetHistogramAttributes(dndetaForward, kRed+1, 20, "ALICE Forward");
  SetHistogramAttributes(normForward, kRed+1, 20, "ALICE Forward", "Normalisation");

  fOutput->Add(fTriggers->Clone());
  fOutput->Add(fSNNString->Clone());
  fOutput->Add(fSysString->Clone());
  fOutput->Add(dndetaForward);
  fOutput->Add(normForward);
  fOutput->Add(Rebin(dndetaForward));

  if (fSumForwardMC) { 
    // Get acceptance normalisation from underflow bins 
    TH1D* normForwardMC   = fSumForwardMC->ProjectionX("normForwardMC", 0, 1, "");
    // Project onto eta axis - _ignoring_underflow_bins_!
    TH1D* dndetaForwardMC = fSumForwardMC->ProjectionX("dndetaForwardMC", 1, -1, "e");
    // Normalize to the acceptance 
    dndetaForwardMC->Divide(normForwardMC);
    // Scale by the vertex efficiency 
    dndetaForwardMC->Scale(vtxEff, "width");

    SetHistogramAttributes(dndetaForwardMC, kRed+3, 21, "ALICE Forward (MC)");
    SetHistogramAttributes(normForwardMC, kRed+3, 21, "ALICE Forward (MC)", 
			   "Normalisation");

    fOutput->Add(dndetaForwardMC);
    fOutput->Add(normForwardMC);
    fOutput->Add(Rebin(dndetaForwardMC));
  }

  if (fSumPrimary) { 
    TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",-1,-1,"e");
    dndetaTruth->Scale(1./nAll, "width");

    SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");

    fOutput->Add(dndetaTruth);
    fOutput->Add(Rebin(dndetaTruth));
  }
  if (fSumCentral) { 
    TH1D* dndetaCentral = fSumCentral->ProjectionX("dndetaCentral",-1,-1,"e");
    dndetaCentral->Scale(1./nGood, "width");

    SetHistogramAttributes(dndetaCentral, kGray+3, 22, "ALICE Central - track(let)s");

    dndetaCentral->GetXaxis()->SetRangeUser(-1,1);
    fOutput->Add(dndetaCentral);
    fOutput->Add(Rebin(dndetaCentral));
  }

  TNamed* trigString = 
    new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
  trigString->SetUniqueID(fTriggerMask);
  fOutput->Add(trigString);

  TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
  vtxAxis->SetName("vtxAxis");
  vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
  fOutput->Add(vtxAxis);

  // If only there was a away to get sqrt{s_NN} and beam type 
  

  PostData(2, fOutput);
}

//________________________________________________________________________
TH1D*
AliForwarddNdetaTask::Rebin(const TH1D* h) const
{
  if (fRebin <= 1) return 0;

  Int_t nBins = h->GetNbinsX();
  if(nBins % fRebin != 0) {
    Warning("Rebin", "Rebin factor %d is not a devisor of current number "
	    "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
    return 0;
  }
    
  // Make a copy 
  TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d", 
					       h->GetName(), fRebin)));
  tmp->Rebin(fRebin);
  tmp->SetDirectory(0);

  // The new number of bins 
  Int_t nBinsNew = nBins / fRebin;
  for(Int_t i = 1;i<= nBinsNew; i++) {
    Double_t content = 0;
    Double_t sumw    = 0;
    Double_t wsum    = 0;
    Int_t    nbins   = 0;
    for(Int_t j = 1; j<=fRebin;j++) {
      Int_t    bin = (i-1)*fRebin + j;
      Double_t c   =  h->GetBinContent(bin);

      if (c <= 0) continue;

      if (fCutEdges) {
	if (h->GetBinContent(bin+1)<=0 || 
	    h->GetBinContent(bin-1)) {
	  Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)", 
		  bin, c, h->GetName(), 
		  bin+1, h->GetBinContent(bin+1), 
		  bin-1, h->GetBinContent(bin-1));
	  continue;
	}	
      }
      Double_t e =  h->GetBinError(bin);
      Double_t w =  1 / (e*e); // 1/c/c
      content    += c;
      sumw       += w;
      wsum       += w * c;
      nbins++;
    }
      
    if(content > 0 && nbins > 0) {
      tmp->SetBinContent(i, wsum / sumw);
      tmp->SetBinError(i,1./TMath::Sqrt(sumw));
    }
  }
  
  return tmp;
}
Commit	Line	Data
b2e7f2d6	1	//====================================================================
	2	#include "AliForwarddNdetaTask.h"
	3	#include <TMath.h>
	4	#include <TH2D.h>
	5	#include <TH1D.h>
	6	#include <THStack.h>
	7	#include <TList.h>
	8	#include <AliAnalysisManager.h>
	9	#include <AliAODEvent.h>
	10	#include <AliAODHandler.h>
	11	#include <AliAODInputHandler.h>
	12	#include "AliForwardUtil.h"
ec82df5b	13	#include <AliAODForwardMult.h>
b2e7f2d6	14
	15	//____________________________________________________________________
	16	AliForwarddNdetaTask::AliForwarddNdetaTask()
	17	: AliAnalysisTaskSE(),
	18	fSumForward(0), // Sum of histograms
	19	fSumForwardMC(0), // Sum of MC histograms (if any)
	20	fSumPrimary(0), // Sum of primary histograms
	21	fSumCentral(0), // Sum of central histograms
	22	fCentral(0), // Cache of central histogram
	23	fPrimary(0), // Cache of primary histogram
	24	fSums(0), // Container of sums
	25	fOutput(0), // Container of outputs
	26	fTriggers(0), // Histogram of triggers
	27	fVtxMin(0), // Minimum v_z
	28	fVtxMax(0), // Maximum v_z
	29	fTriggerMask(0), // Trigger mask
	30	fRebin(0), // Rebinning factor
	31	fCutEdges(false),
	32	fSNNString(0),
	33	fSysString(0)
	34	{}
	35
	36	//____________________________________________________________________
	37	AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
	38	: AliAnalysisTaskSE("Forward"),
	39	fSumForward(0), // Sum of histograms
	40	fSumForwardMC(0), // Sum of MC histograms (if any)
	41	fSumPrimary(0), // Sum of primary histograms
	42	fSumCentral(0), // Sum of central histograms
	43	fCentral(0), // Cache of central histogram
	44	fPrimary(0), // Cache of primary histogram
	45	fSums(0), // Container of sums
	46	fOutput(0), // Container of outputs
	47	fTriggers(0), // Histogram of triggers
	48	fVtxMin(-10), // Minimum v_z
	49	fVtxMax(10), // Maximum v_z
	50	fTriggerMask(AliAODForwardMult::kInel),
	51	fRebin(5), // Rebinning factor
	52	fCutEdges(false),
	53	fSNNString(0),
	54	fSysString(0)
	55	{
	56	// Output slot #1 writes into a TH1 container
	57	DefineOutput(1, TList::Class());
	58	DefineOutput(2, TList::Class());
	59	}
	60
	61	//____________________________________________________________________
	62	AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
	63	: AliAnalysisTaskSE(o),
ec82df5b	64	fSumForward(o.fSumForward), // Sum of histograms
	65	fSumForwardMC(o.fSumForwardMC),// Sum of MC histograms (if any)
	66	fSumPrimary(o.fSumPrimary), // Sum of primary histograms
	67	fSumCentral(o.fSumCentral), // Sum of central histograms
	68	fCentral(o.fCentral), // Cache of central histogram
	69	fPrimary(o.fPrimary), // Cache of primary histogram
	70	fSums(o.fSums), // Container of sums
	71	fOutput(o.fOutput), // Container of outputs
	72	fTriggers(o.fTriggers), // Histogram of triggers
	73	fVtxMin(o.fVtxMin), // Minimum v_z
	74	fVtxMax(o.fVtxMax), // Maximum v_z
	75	fTriggerMask(o.fTriggerMask), // Trigger mask
	76	fRebin(o.fRebin), // Rebinning factor
	77	fCutEdges(o.fCutEdges), // Whether to cut edges when rebinning
	78	fSNNString(o.fSNNString), //
	79	fSysString(o.fSysString) //
b2e7f2d6	80	{}
	81
	82	//____________________________________________________________________
	83	AliForwarddNdetaTask::~AliForwarddNdetaTask()
	84	{
	85	if (fSums) {
	86	fSums->Delete();
	87	delete fSums;
	88	fSums = 0;
	89	}
	90	if (fOutputs) {
	91	fOutputs->Delete();
	92	delete fOutputs;
	93	fOutputs = 0;
	94	}
	95	if (fCentral) delete fCentral;
	96	if (fPrimary) delete fPrimary;
	97	}
	98
	99	//________________________________________________________________________
	100	void
	101	AliForwarddNdetaTask::SetTriggerMask(const char* mask)
	102	{
	103	UShort_t trgMask = 0;
	104	TString trgs(mask);
	105	trgs.ToUpper();
	106	TObjString* trg;
	107	TIter next(trgs.Tokenize(" ,\|"));
	108	while ((trg = static_cast<TObjString*>(next()))) {
	109	TString s(trg->GetString());
	110	if (s.IsNull()) continue;
	111	if (s.CompareTo("INEL") == 0) trgMask = AliAODForwardMult::kInel;
	112	else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
	113	else if (s.CompareTo("NSD") == 0) trgMask = AliAODForwardMult::kNSD;
	114	else
	115	Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
	116	}
	117	if (trgMask == 0) trgMask = 1;
	118	SetTriggerMask(trgMask);
	119	}
	120
	121	//________________________________________________________________________
	122	void
	123	AliForwarddNdetaTask::UserCreateOutputObjects()
	124	{
	125	// Create histograms
	126	// Called once (on the worker node)
	127
	128	fOutput = new TList;
	129	fOutput->SetName(Form("%s_result", GetName()));
	130	fOutput->SetOwner();
	131
	132	fSums = new TList;
	133	fSums->SetName(Form("%s_sums", GetName()));
	134	fSums->SetOwner();
	135
	136
	137	fTriggers = new TH1D("triggers", "Number of triggers",
	138	kAccepted, 1, kAccepted);
	139	fTriggers->SetYTitle("# of events");
	140	fTriggers->GetXaxis()->SetBinLabel(kAll, "All events");
	141	fTriggers->GetXaxis()->SetBinLabel(kB, "w/B trigger");
	142	fTriggers->GetXaxis()->SetBinLabel(kA, "w/A trigger");
	143	fTriggers->GetXaxis()->SetBinLabel(kC, "w/C trigger");
144	fTriggers->GetXaxis()->SetBinLabel(kE, "w/E trigger");
145	fTriggers->GetXaxis()->SetBinLabel(kMB, "w/Collision trigger");
146	fTriggers->GetXaxis()->SetBinLabel(kWithVertex, "w/Vertex");
147	fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
148	fTriggers->GetXaxis()->SetBinLabel(kAccepted, "Accepted by cut");
149	fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
150	fTriggers->SetFillColor(kRed+1);
151	fTriggers->SetFillStyle(3001);
152	fTriggers->SetStats(0);
153	fSums->Add(fTriggers);
154
155	fSNNString = new TNamed("sNN", "");
156	fSysString = new TNamed("sys", "");
157	fSums->Add(fSNNString);
158	fSums->Add(fSysString);
159
160	// Check that we have an AOD input handler
161	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
162	AliAODInputHandler* ah =
163	dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
164	if (!ah) AliFatal("No AOD input handler set in analysis manager");
165
166	// Post data for ALL output slots >0 here, to get at least an empty histogram
167	PostData(1, fSums);
168	}
169
170	//____________________________________________________________________
171	TH2D*
172	AliForwarddNdetaTask::CloneHist(TH2D* in, const char* name)
173	{
174	if (!in) return 0;
175	TH2D* ret = static_cast<TH2D*>(in->Clone(name));
176	ret->SetDirectory(0);
177	ret->Sumw2();
178	fSums->Add(ret);
179
180	return ret;
181	}
182
183	//____________________________________________________________________
184	void
185	AliForwarddNdetaTask::UserExec(Option_t *)
186	{
187	// Main loop
188	AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
189	if (!aod) {
190	AliError("Cannot get the AOD event");
191	return;
192	}
193
194	// Get objects from the event structure
195	TObject* oForward = aod->FindListObject("Forward");
196	TObject* oForwardMC = aod->FindListObject("ForwardMC");
197	TObject* oPrimary = aod->FindListObject("primary");
198	TObject* oCentral = aod->FindListObject("Central");
199
200	// We should have a forward object at least
201	if (!oForward) {
202	AliWarning("No Forward object found AOD");
203	return;
204	}
205
206	// Cast to good types
207	AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(oForward);
208	AliAODForwardMult* forwardMC = static_cast<AliAODForwardMult*>(oForwardMC);
209	TH2D* primary = static_cast<TH2D*>(oPrimary);
210	TH2D* central = static_cast<TH2D*>(oCentral);
211
212	static bool first = true;
213	if (first) {
214	UShort_t sNN = forward->GetSNN();
215	UShort_t sys = forward->GetSystem();
216	fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
217	fSNNString->SetUniqueID(sNN);
218	fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
219	fSysString->SetUniqueID(sys);
220	first = false;
221	}
222
223	// Create our sum histograms
224	if (!fSumForward) fSumForward = CloneHist(&forward->GetHistogram(),"forward");
225	if (!fSumForwardMC && forwardMC)
226	fSumForwardMC = CloneHist(&forwardMC->GetHistogram(),"forwardMC");
227	if (!fSumPrimary && primary) fSumPrimary = CloneHist(primary, "primary");
228	if (!fSumCentral && central) fSumCentral = CloneHist(central, "central");
229
230	// Add contribtion from MC
231	if (primary) fSumPrimary->Add(primary);
232
233	// Count event
234	fTriggers->AddBinContent(kAll);
235	if (forward->IsTriggerBits(AliAODForwardMult::kB))
236	fTriggers->AddBinContent(kB);
237	if (forward->IsTriggerBits(AliAODForwardMult::kA))
238	fTriggers->AddBinContent(kA);
239	if (forward->IsTriggerBits(AliAODForwardMult::kC))
240	fTriggers->AddBinContent(kC);
241	if (forward->IsTriggerBits(AliAODForwardMult::kE))
242	fTriggers->AddBinContent(kE);
243	if (forward->IsTriggerBits(AliAODForwardMult::kInel))
244	fTriggers->AddBinContent(kMB);
245
246	// Check if we have an event of interest.
247	if (!forward->IsTriggerBits(fTriggerMask)) return;
248	fTriggers->AddBinContent(kWithTrigger);
249
250	// Check that we have a valid vertex
251	if (!forward->HasIpZ()) return;
252	fTriggers->AddBinContent(kWithVertex);
253
254	// Check that vertex is within cuts
255	if (!forward->InRange(fVtxMin, fVtxMax)) return;
256	fTriggers->AddBinContent(kAccepted);
257
258	// Add contribution
259	fSumForward->Add(&(forward->GetHistogram()));
260	if (fSumForwardMC) fSumForwardMC->Add(&(forwardMC->GetHistogram()));
261	if (fSumPrimary) fSumPrimary->Add(primary);
262	if (fSumCentral) fSumCentral->Add(central);
263
264	PostData(1, fSums);
265	}
266
267	//________________________________________________________________________
268	void
269	AliForwarddNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
270	const char* title, const char* ytitle)
271	{
272	h->SetTitle(title);
273	h->SetMarkerColor(colour);
274	h->SetMarkerStyle(marker);
275	h->SetMarkerSize(1);
276	h->SetFillStyle(0);
277	h->SetYTitle(ytitle);
278	h->GetXaxis()->SetTitleFont(132);
279	h->GetXaxis()->SetLabelFont(132);
280	h->GetXaxis()->SetNdivisions(10);
281	h->GetYaxis()->SetTitleFont(132);
282	h->GetYaxis()->SetLabelFont(132);
283	h->GetYaxis()->SetNdivisions(10);
284	h->GetYaxis()->SetDecimals();
285	h->SetStats(0);
286	}
287
288	//________________________________________________________________________
289	void
290	AliForwarddNdetaTask::Terminate(Option_t *)
291	{
292	// Draw result to screen, or perform fitting, normalizations Called
293	// once at the end of the query
294
295	fSums = dynamic_cast<TList*> (GetOutputData(1));
8e2bb72a	296	if(!fSums) {
b2e7f2d6	297	AliError("Could not retrieve TList fSums");
	298	return;
	299	}
	300
	301	if (!fOutput) {
	302	fOutput = new TList;
	303	fOutput->SetName(Form("%s_result", GetName()));
	304	fOutput->SetOwner();
	305	}
	306
	307	fSumForward = static_cast<TH2D*>(fSums->FindObject("forward"));
	308	fSumForwardMC = static_cast<TH2D*>(fSums->FindObject("forwardMC"));
	309	fSumPrimary = static_cast<TH2D*>(fSums->FindObject("primary"));
	310	fSumCentral = static_cast<TH2D*>(fSums->FindObject("central"));
	311	fTriggers = static_cast<TH1D*>(fSums->FindObject("triggers"));
	312	fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
	313	fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
	314
	315	if (!fTriggers) {
	316	AliError("Couldn't find histogram 'triggers' in list");
	317	return;
	318	}
	319	if (!fSumForward) {
	320	AliError("Couldn't find histogram 'forward' in list");
	321	return;
	322	}
	323
	324	Int_t nAll = Int_t(fTriggers->GetBinContent(kAll));
	325	Int_t nB = Int_t(fTriggers->GetBinContent(kB));
	326	Int_t nA = Int_t(fTriggers->GetBinContent(kA));
	327	Int_t nC = Int_t(fTriggers->GetBinContent(kC));
	328	Int_t nE = Int_t(fTriggers->GetBinContent(kE));
	329	Int_t nMB = Int_t(fTriggers->GetBinContent(kMB));
	330	Int_t nTriggered = Int_t(fTriggers->GetBinContent(kWithTrigger));
	331	Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
	332	Int_t nAccepted = Int_t(fTriggers->GetBinContent(kAccepted));
	333	Int_t nGood = nB - nA - nC + 2 * nE;
	334	Double_t vtxEff = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
	335	AliInfo(Form("Total of %9d events\n"
	336	" of these %9d are minimum bias\n"
	337	" of these %9d has a %s trigger\n"
	338	" of these %9d has a vertex\n"
	339	" of these %9d were in [%+4.1f,%+4.1f]cm\n"
	340	" Triggers by type:\n"
	341	" B = %9d\n"
	342	" A\|C = %9d (%9d+%-9d)\n"
	343	" E = %9d\n"
	344	" Implies %9d good triggers\n"
	345	" Vertex efficiency: %f",
	346	nAll, nMB, nTriggered,
	347	AliAODForwardMult::GetTriggerString(fTriggerMask),
	348	nWithVertex, nAccepted,
	349	fVtxMin, fVtxMax,
	350	nB, nA+nC, nA, nC, nE, nGood, vtxEff));
	351
	352	// Get acceptance normalisation from underflow bins
	353	TH1D* normForward = fSumForward->ProjectionX("normForward", 0, 1, "");
	354	// Project onto eta axis - _ignoring_underflow_bins_!
	355	TH1D* dndetaForward = fSumForward->ProjectionX("dndetaForward", 1, -1, "e");
	356	// Normalize to the acceptance
	357	dndetaForward->Divide(normForward);
	358	// Scale by the vertex efficiency
	359	dndetaForward->Scale(vtxEff, "width");
	360
361	SetHistogramAttributes(dndetaForward, kRed+1, 20, "ALICE Forward");
362	SetHistogramAttributes(normForward, kRed+1, 20, "ALICE Forward", "Normalisation");
363
364	fOutput->Add(fTriggers->Clone());
365	fOutput->Add(fSNNString->Clone());
366	fOutput->Add(fSysString->Clone());
367	fOutput->Add(dndetaForward);
368	fOutput->Add(normForward);
369	fOutput->Add(Rebin(dndetaForward));
370
371	if (fSumForwardMC) {
372	// Get acceptance normalisation from underflow bins
373	TH1D* normForwardMC = fSumForwardMC->ProjectionX("normForwardMC", 0, 1, "");
374	// Project onto eta axis - _ignoring_underflow_bins_!
375	TH1D* dndetaForwardMC = fSumForwardMC->ProjectionX("dndetaForwardMC", 1, -1, "e");
376	// Normalize to the acceptance
377	dndetaForwardMC->Divide(normForwardMC);
378	// Scale by the vertex efficiency
379	dndetaForwardMC->Scale(vtxEff, "width");
380
381	SetHistogramAttributes(dndetaForwardMC, kRed+3, 21, "ALICE Forward (MC)");
382	SetHistogramAttributes(normForwardMC, kRed+3, 21, "ALICE Forward (MC)",
383	"Normalisation");
384
385	fOutput->Add(dndetaForwardMC);
386	fOutput->Add(normForwardMC);
387	fOutput->Add(Rebin(dndetaForwardMC));
388	}
389
390	if (fSumPrimary) {
391	TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",-1,-1,"e");
392	dndetaTruth->Scale(1./nAll, "width");
393
394	SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");
395
396	fOutput->Add(dndetaTruth);
397	fOutput->Add(Rebin(dndetaTruth));
398	}
399	if (fSumCentral) {
400	TH1D* dndetaCentral = fSumCentral->ProjectionX("dndetaCentral",-1,-1,"e");
401	dndetaCentral->Scale(1./nGood, "width");
402
403	SetHistogramAttributes(dndetaCentral, kGray+3, 22, "ALICE Central - track(let)s");
404
405	dndetaCentral->GetXaxis()->SetRangeUser(-1,1);
406	fOutput->Add(dndetaCentral);
407	fOutput->Add(Rebin(dndetaCentral));
408	}
409
410	TNamed* trigString =
411	new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
412	trigString->SetUniqueID(fTriggerMask);
413	fOutput->Add(trigString);
414
415	TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
416	vtxAxis->SetName("vtxAxis");
417	vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
418	fOutput->Add(vtxAxis);
419
420	// If only there was a away to get sqrt{s_NN} and beam type
421
422
423	PostData(2, fOutput);
424	}
425
426	//________________________________________________________________________
427	TH1D*
428	AliForwarddNdetaTask::Rebin(const TH1D* h) const
429	{
430	if (fRebin <= 1) return 0;
431
432	Int_t nBins = h->GetNbinsX();
433	if(nBins % fRebin != 0) {
434	Warning("Rebin", "Rebin factor %d is not a devisor of current number "
435	"of bins %d in the histogram %s", fRebin, nBins, h->GetName());
436	return 0;
437	}
438
439	// Make a copy
440	TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
441	h->GetName(), fRebin)));
442	tmp->Rebin(fRebin);
443	tmp->SetDirectory(0);
444
445	// The new number of bins
446	Int_t nBinsNew = nBins / fRebin;
447	for(Int_t i = 1;i<= nBinsNew; i++) {
448	Double_t content = 0;
449	Double_t sumw = 0;
450	Double_t wsum = 0;
451	Int_t nbins = 0;
452	for(Int_t j = 1; j<=fRebin;j++) {
453	Int_t bin = (i-1)*fRebin + j;
454	Double_t c = h->GetBinContent(bin);
455
456	if (c <= 0) continue;
457
458	if (fCutEdges) {
459	if (h->GetBinContent(bin+1)<=0 \|\|
460	h->GetBinContent(bin-1)) {
461	Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
462	bin, c, h->GetName(),
463	bin+1, h->GetBinContent(bin+1),
464	bin-1, h->GetBinContent(bin-1));
465	continue;
466	}
467	}
468	Double_t e = h->GetBinError(bin);
469	Double_t w = 1 / (e*e); // 1/c/c
470	content += c;
471	sumw += w;
472	wsum += w * c;
473	nbins++;
474	}
475
476	if(content > 0 && nbins > 0) {
477	tmp->SetBinContent(i, wsum / sumw);
478	tmp->SetBinError(i,1./TMath::Sqrt(sumw));
479	}
480	}
481
482	return tmp;
483	}