[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 
}

//____________________________________________________________________
AliMCTruthdNdetaTask::AliMCTruthdNdetaTask(const AliMCTruthdNdetaTask& o)
  : AliBasedNdetaTask(o), 
    fHasData(o.fHasData)
{
  // 
  // Copy constructor
  // 
}

//____________________________________________________________________
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;
  TObject* obj = aod->FindListObject("primary");
  // We should have a forward object at least 
  if (!obj) {
    fHasData = false;
    return 0;
  }
  TH2D* ret = static_cast<TH2D*>(obj);
  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;
}

//________________________________________________________________________
void 
AliMCTruthdNdetaTask::Terminate(Option_t *option) 
{
  // 
  // 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;
  }
  AliBasedNdetaTask::Terminate(option);

  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 (fCentAxis && 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) fOutput->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) fOutput->Add(truthRebin);  
}

//========================================================================
void
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;

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

  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);
}

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