[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliMCTruthdNdetaTask.cxx

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

//____________________________________________________________________
AliMCTruthdNdetaTask::AliMCTruthdNdetaTask()
  : AliBasedNdetaTask(),
    fHasData(true)
{
  //
  // Constructor 
  // 
}

//____________________________________________________________________
AliMCTruthdNdetaTask::AliMCTruthdNdetaTask(const char* /* name */)
  : AliBasedNdetaTask("MCTruth"), 
    fHasData(true)
{
  // 
  // Constructor
  // 
  // Paramters
  //   name    Name of task 
}


//____________________________________________________________________
AliBasedNdetaTask::CentralityBin*
AliMCTruthdNdetaTask::MakeCentralityBin(const char* name, Short_t l, Short_t h) 
  const 
{
  // 
  // Make a new centrality bin
  // 
  // Parameters:
  //    name   Histogram names
  //    l      Lower cut
  //    h      Upper cut
  // 
  // Return:
  //    Newly allocated object (of our type)
  //
  return new AliMCTruthdNdetaTask::CentralityBin(name, l, h);
}

//____________________________________________________________________
TH2D*
AliMCTruthdNdetaTask::GetHistogram(const AliAODEvent& aod, Bool_t mc)
{
  // 
  // Retrieve the histogram 
  // 
  // Parameters:
  //    aod AOD event 
  //    mc  Whether to get the MC histogram or not
  // 
  // Return:
  //    Retrieved histogram or null
  //
  if (!fHasData) return 0;
  if (mc) return 0;

  TH2D* ret = GetPrimary(aod);
  if (!ret) {
    fHasData = false;
    return 0;
  }
  Int_t nY  = ret->GetNbinsY();
  // Need to fill under-/overflow bin with 1's 
  for (Int_t i = 1; i <= ret->GetNbinsX(); i++)  {
    ret->SetBinContent(i, 0,    1);
    ret->SetBinContent(i, nY+1, 1);
  }
  return ret;
}

//________________________________________________________________________
Bool_t
AliMCTruthdNdetaTask::Finalize() 
{
  // 
  // Called at end of event processing.. 
  //
  // This is called once in the master 
  // 
  // Parameters:
  //    option Not used 
  if (!fHasData) {
    AliInfo("The MC truth dN/deta task didn't get any data");
    return false;
  }
  AliBasedNdetaTask::Finalize();

  THStack* truth      = new THStack("dndetaTruth", "dN/d#eta MC Truth");
  THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin), 
				    "dN/d#eta MC Truth");

  TIter next(fListOfCentralities);
  CentralityBin* bin = 0;
  while ((bin = static_cast<CentralityBin*>(next()))) {
    if (HasCentrality() && bin->IsAllBin()) continue;

    TList* results = bin->GetResults();
    if (!results) continue; 

    TH1* dndeta      = static_cast<TH1*>(results->FindObject("dndetaTruth"));
    TH1* dndetaRebin = 
      static_cast<TH1*>(results->FindObject(Form("dndetaTruth_rebin%02d",
						fRebin)));
    if (dndeta)      truth->Add(dndeta);
    if (dndetaRebin) truthRebin->Add(dndetaRebin);
  }
  // If available output rebinned stack 
  if (!truth->GetHists() || 
      truth->GetHists()->GetEntries() <= 0) {
    AliWarning("No MC truth histograms found");
    delete truth;
    truth = 0;
  }
  if (truth) fResults->Add(truth);

  // If available output rebinned stack 
  if (!truthRebin->GetHists() || 
      truthRebin->GetHists()->GetEntries() <= 0) {
    AliWarning("No rebinned MC truth histograms found");
    delete truthRebin;
    truthRebin = 0;
  }
  if (truthRebin) fResults->Add(truthRebin);  
  return true;
}

//========================================================================
Bool_t
AliMCTruthdNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* 
						  forward, 
						  Int_t triggerMask,
						  Bool_t isZero,
						  Double_t vzMin, 
						  Double_t vzMax, 
						  const TH2D* primary,
						  const TH2D*)
{ 
  // Check the centrality class unless this is the 'all' bin 
  if (!primary) return false;

  if (!IsAllBin()) { 
    Double_t centrality = forward->GetCentrality();
    if (centrality < fLow || centrality >= fHigh) return false;
  }

  if (!fSum) CreateSums(primary, 0);
  if (!fSumTruth) { 
    fSumTruth = static_cast<TH2D*>(primary->Clone("truth"));
    fSumTruth->SetDirectory(0);
    fSumTruth->Reset();
    fSums->Add(fSumTruth);
  }

  // translate real trigger mask to MC trigger mask
  Int_t mask = AliAODForwardMult::kB;
  if (triggerMask == AliAODForwardMult::kNSD) {
    mask ^= AliAODForwardMult::kNSD;
    mask =  AliAODForwardMult::kMCNSD;
  }

  // Now use our normal check, but with the new mask, except ignore vertex
  if (forward->CheckEvent(mask, -10000, -10000, 0, 0, 0)) {
    fSumTruth->Add(primary);

    // Store event count in left-most under- underflow bin 
    Int_t cnt = Int_t(fSumTruth->GetBinContent(0,0));
    fSumTruth->SetBinContent(0,0, ++cnt);
  }

  // Now use our normal check with the full trigger mask and vertex
  if (CheckEvent(forward, triggerMask, vzMin, vzMax)) 
    fSum->Add(primary, isZero);
  return true;
}

//________________________________________________________________________
void
AliMCTruthdNdetaTask::CentralityBin::End(TList*      sums, 
					 TList*      results,
					 UShort_t    scheme,
					 const TH2F* shapeCorr, 
					 Double_t    trigEff,
					 Double_t    trigEff0,
					 Bool_t      symmetrice,
					 Int_t       rebin, 
					 Bool_t      rootProj,
					 Bool_t      corrEmpty, 
					 Bool_t      cutEdges,
					 Int_t       triggerMask,
					 Int_t       marker,
					 Int_t       color,
					 TList*      mclist,
					 TList*      truthlist)
{
#if 0
  AliInfo(Form("At end with sums=%p, results=%p, scheme=%d, "
	       "shapeCorr=%p, trigEff=%f, symmetrice=%d, rebin=%d, "
	       "rootProj=%d, corrEmpty=%d, cutEdges=%d, triggerMask=0x%08x, "
	       "marker=%d (%d)", 
	       sums, results, scheme, shapeCorr, trigEff, trigEff0, symmetrice, 
	       rebin, rootProj, corrEmpty, cutEdges, triggerMask, marker,
	       GetMarkerStyle(kStar)));
#endif

  AliBasedNdetaTask::CentralityBin::End(sums, results, scheme, 
					shapeCorr, trigEff, trigEff0,
					symmetrice, rebin, 
					rootProj, corrEmpty, cutEdges,
					triggerMask, marker, color, mclist, 
					truthlist);

  fSumTruth     = static_cast<TH2D*>(fSums->FindObject("truth"));
  

  if (fSumTruth) { 
    Int_t n0 = Int_t(fSumTruth->GetBinContent(0,0));
    Int_t n  = (triggerMask == AliAODForwardMult::kNSD ? 
		Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) : 
		Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
    DMSG(fDebug,0,"Normalising MC truth to %d (%d additions)", n, n0);
    
    TH1D* dndetaTruth = fSumTruth->ProjectionX("dndetaTruth",1,
					       fSumTruth->GetNbinsY(),"e");
    dndetaTruth->SetDirectory(0);
    dndetaTruth->Scale(1./n, "width");
    
    SetHistogramAttributes(dndetaTruth, GetColor(color)+1, 
			   GetMarkerStyle(kCross), 
			   "Monte-Carlo truth");

    fOutput->Add(dndetaTruth);
    fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
  }
  TH1* dndeta          =                       GetResult(0,     false, "");
  TH1* dndetaSym       = symmetrice          ? GetResult(0,     true,  "") : 0;
  TH1* dndeta_rebin    = rebin               ? GetResult(rebin, false, "") : 0;
  TH1* dndetaSym_rebin = symmetrice && rebin ? GetResult(rebin, true,  "") : 0;
  if (dndeta)    
    dndeta->SetTitle("Monte-Carlo truth (selected)");
  if (dndetaSym) 
    dndetaSym->SetTitle("Monte-Carlo truth (selected,mirrored)");
  if (dndeta_rebin)    
    dndeta_rebin->SetTitle("Monte-Carlo truth (selected)");
  if (dndetaSym_rebin) 
    dndetaSym_rebin->SetTitle("Monte-Carlo truth (selected,mirrored)");

}

//________________________________________________________________________
//
// EOF
//
Commit	Line	Data
9ecab72f	1	//====================================================================
	2	#include "AliMCTruthdNdetaTask.h"
	3	#include <TMath.h>
	4	#include <TH2D.h>
	5	#include <TH1D.h>
	6	#include <THStack.h>
	7	#include <TList.h>
	8	#include <TFile.h>
	9	#include <AliAnalysisManager.h>
	10	#include <AliAODEvent.h>
	11	#include <AliAODHandler.h>
	12	#include <AliAODInputHandler.h>
	13	#include "AliForwardUtil.h"
	14	#include "AliAODForwardMult.h"
	15
	16	//____________________________________________________________________
	17	AliMCTruthdNdetaTask::AliMCTruthdNdetaTask()
786db950	18	: AliBasedNdetaTask(),
786db950	19	fHasData(true)
9ecab72f	20	{
	21	//
	22	// Constructor
	23	//
	24	}
	25
	26	//____________________________________________________________________
	27	AliMCTruthdNdetaTask::AliMCTruthdNdetaTask(const char* /* name */)
786db950	28	: AliBasedNdetaTask("MCTruth"),
786db950	29	fHasData(true)
9ecab72f	30	{
	31	//
	32	// Constructor
	33	//
	34	// Paramters
	35	// name Name of task
	36	}
	37
9ecab72f	38
	39	//____________________________________________________________________
	40	AliBasedNdetaTask::CentralityBin*
	41	AliMCTruthdNdetaTask::MakeCentralityBin(const char* name, Short_t l, Short_t h)
	42	const
	43	{
	44	//
	45	// Make a new centrality bin
	46	//
	47	// Parameters:
	48	// name Histogram names
	49	// l Lower cut
	50	// h Upper cut
	51	//
	52	// Return:
	53	// Newly allocated object (of our type)
	54	//
	55	return new AliMCTruthdNdetaTask::CentralityBin(name, l, h);
	56	}
	57
	58	//____________________________________________________________________
	59	TH2D*
c8b1a7db	60	AliMCTruthdNdetaTask::GetHistogram(const AliAODEvent& aod, Bool_t mc)
9ecab72f	61	{
	62	//
	63	// Retrieve the histogram
	64	//
	65	// Parameters:
	66	// aod AOD event
	67	// mc Whether to get the MC histogram or not
	68	//
	69	// Return:
	70	// Retrieved histogram or null
	71	//
786db950	72	if (!fHasData) return 0;
9ecab72f	73	if (mc) return 0;
c8b1a7db	74
	75	TH2D* ret = GetPrimary(aod);
	76	if (!ret) {
786db950	77	fHasData = false;
	78	return 0;
	79	}
808d97d9	80	Int_t nY = ret->GetNbinsY();
	81	// Need to fill under-/overflow bin with 1's
	82	for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
	83	ret->SetBinContent(i, 0, 1);
	84	ret->SetBinContent(i, nY+1, 1);
	85	}
9ecab72f	86	return ret;
	87	}
	88
	89	//________________________________________________________________________
c8b1a7db	90	Bool_t
c8b1a7db	91	AliMCTruthdNdetaTask::Finalize()
9ecab72f	92	{
	93	//
	94	// Called at end of event processing..
	95	//
	96	// This is called once in the master
	97	//
	98	// Parameters:
	99	// option Not used
786db950	100	if (!fHasData) {
786db950	101	AliInfo("The MC truth dN/deta task didn't get any data");
c8b1a7db	102	return false;
786db950	103	}
c8b1a7db	104	AliBasedNdetaTask::Finalize();
9ecab72f	105
	106	THStack* truth = new THStack("dndetaTruth", "dN/d#eta MC Truth");
	107	THStack* truthRebin = new THStack(Form("dndetaTruth_rebin%02d", fRebin),
	108	"dN/d#eta MC Truth");
	109
	110	TIter next(fListOfCentralities);
	111	CentralityBin* bin = 0;
	112	while ((bin = static_cast<CentralityBin*>(next()))) {
c8b1a7db	113	if (HasCentrality() && bin->IsAllBin()) continue;
9ecab72f	114
	115	TList* results = bin->GetResults();
	116	if (!results) continue;
	117
	118	TH1* dndeta = static_cast<TH1*>(results->FindObject("dndetaTruth"));
	119	TH1* dndetaRebin =
	120	static_cast<TH1*>(results->FindObject(Form("dndetaTruth_rebin%02d",
	121	fRebin)));
	122	if (dndeta) truth->Add(dndeta);
	123	if (dndetaRebin) truthRebin->Add(dndetaRebin);
	124	}
	125	// If available output rebinned stack
	126	if (!truth->GetHists() \|\|
	127	truth->GetHists()->GetEntries() <= 0) {
	128	AliWarning("No MC truth histograms found");
	129	delete truth;
	130	truth = 0;
	131	}
c8b1a7db	132	if (truth) fResults->Add(truth);
9ecab72f	133
	134	// If available output rebinned stack
	135	if (!truthRebin->GetHists() \|\|
	136	truthRebin->GetHists()->GetEntries() <= 0) {
	137	AliWarning("No rebinned MC truth histograms found");
	138	delete truthRebin;
	139	truthRebin = 0;
	140	}
c8b1a7db	141	if (truthRebin) fResults->Add(truthRebin);
c8b1a7db	142	return true;
9ecab72f	143	}
	144
	145	//========================================================================
bfab35d9	146	Bool_t
9ecab72f	147	AliMCTruthdNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult*
	148	forward,
	149	Int_t triggerMask,
	150	Bool_t isZero,
	151	Double_t vzMin,
	152	Double_t vzMax,
	153	const TH2D* primary,
	154	const TH2D*)
	155	{
	156	// Check the centrality class unless this is the 'all' bin
bfab35d9	157	if (!primary) return false;
9ecab72f	158
	159	if (!IsAllBin()) {
	160	Double_t centrality = forward->GetCentrality();
bfab35d9	161	if (centrality < fLow \|\| centrality >= fHigh) return false;
9ecab72f	162	}
	163
	164	if (!fSum) CreateSums(primary, 0);
	165	if (!fSumTruth) {
	166	fSumTruth = static_cast<TH2D*>(primary->Clone("truth"));
	167	fSumTruth->SetDirectory(0);
	168	fSumTruth->Reset();
	169	fSums->Add(fSumTruth);
	170	}
	171
	172	// translate real trigger mask to MC trigger mask
	173	Int_t mask = AliAODForwardMult::kB;
	174	if (triggerMask == AliAODForwardMult::kNSD) {
	175	mask ^= AliAODForwardMult::kNSD;
	176	mask = AliAODForwardMult::kMCNSD;
	177	}
	178
	179	// Now use our normal check, but with the new mask, except ignore vertex
	180	if (forward->CheckEvent(mask, -10000, -10000, 0, 0, 0)) {
	181	fSumTruth->Add(primary);
	182
	183	// Store event count in left-most under- underflow bin
	184	Int_t cnt = Int_t(fSumTruth->GetBinContent(0,0));
	185	fSumTruth->SetBinContent(0,0, ++cnt);
	186	}
	187
	188	// Now use our normal check with the full trigger mask and vertex
	189	if (CheckEvent(forward, triggerMask, vzMin, vzMax))
	190	fSum->Add(primary, isZero);
bfab35d9	191	return true;
9ecab72f	192	}
	193
	194	//________________________________________________________________________
	195	void
	196	AliMCTruthdNdetaTask::CentralityBin::End(TList* sums,
	197	TList* results,
	198	UShort_t scheme,
f67d699c	199	const TH2F* shapeCorr,
9ecab72f	200	Double_t trigEff,
66cf95f2	201	Double_t trigEff0,
9ecab72f	202	Bool_t symmetrice,
	203	Int_t rebin,
	204	Bool_t rootProj,
	205	Bool_t corrEmpty,
	206	Bool_t cutEdges,
	207	Int_t triggerMask,
	208	Int_t marker,
c89b9ac1	209	Int_t color,
	210	TList* mclist,
	211	TList* truthlist)
9ecab72f	212	{
	213	#if 0
	214	AliInfo(Form("At end with sums=%p, results=%p, scheme=%d, "
	215	"shapeCorr=%p, trigEff=%f, symmetrice=%d, rebin=%d, "
	216	"rootProj=%d, corrEmpty=%d, cutEdges=%d, triggerMask=0x%08x, "
	217	"marker=%d (%d)",
66cf95f2	218	sums, results, scheme, shapeCorr, trigEff, trigEff0, symmetrice,
9ecab72f	219	rebin, rootProj, corrEmpty, cutEdges, triggerMask, marker,
	220	GetMarkerStyle(kStar)));
	221	#endif
	222
	223	AliBasedNdetaTask::CentralityBin::End(sums, results, scheme,
66cf95f2	224	shapeCorr, trigEff, trigEff0,
9ecab72f	225	symmetrice, rebin,
9ecab72f	226	rootProj, corrEmpty, cutEdges,
f7cfc454	227	triggerMask, marker, color, mclist,
f7cfc454	228	truthlist);
9ecab72f	229
	230	fSumTruth = static_cast<TH2D*>(fSums->FindObject("truth"));
	231
	232
	233	if (fSumTruth) {
f7cfc454	234	Int_t n0 = Int_t(fSumTruth->GetBinContent(0,0));
9ecab72f	235	Int_t n = (triggerMask == AliAODForwardMult::kNSD ?
	236	Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinMCNSD)) :
	237	Int_t(fTriggers->GetBinContent(AliAODForwardMult::kBinAll)));
c8b1a7db	238	DMSG(fDebug,0,"Normalising MC truth to %d (%d additions)", n, n0);
9ecab72f	239
	240	TH1D* dndetaTruth = fSumTruth->ProjectionX("dndetaTruth",1,
	241	fSumTruth->GetNbinsY(),"e");
	242	dndetaTruth->SetDirectory(0);
	243	dndetaTruth->Scale(1./n, "width");
	244
	245	SetHistogramAttributes(dndetaTruth, GetColor(color)+1,
	246	GetMarkerStyle(kCross),
	247	"Monte-Carlo truth");
	248
	249	fOutput->Add(dndetaTruth);
	250	fOutput->Add(Rebin(dndetaTruth, rebin, cutEdges));
	251	}
c8b1a7db	252	TH1* dndeta = GetResult(0, false, "");
	253	TH1* dndetaSym = symmetrice ? GetResult(0, true, "") : 0;
	254	TH1* dndeta_rebin = rebin ? GetResult(rebin, false, "") : 0;
	255	TH1* dndetaSym_rebin = symmetrice && rebin ? GetResult(rebin, true, "") : 0;
9ecab72f	256	if (dndeta)
	257	dndeta->SetTitle("Monte-Carlo truth (selected)");
	258	if (dndetaSym)
	259	dndetaSym->SetTitle("Monte-Carlo truth (selected,mirrored)");
	260	if (dndeta_rebin)
	261	dndeta_rebin->SetTitle("Monte-Carlo truth (selected)");
	262	if (dndetaSym_rebin)
	263	dndetaSym_rebin->SetTitle("Monte-Carlo truth (selected,mirrored)");
	264
	265	}
	266
	267	//________________________________________________________________________
	268	//
	269	// EOF
	270	//