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

//====================================================================
// 
// Base class for classes that calculate the multiplicity in the
// forward regions event-by-event
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODForwardMult 
// 
// Histograms 
//   
// Corrections used 
#include "AliForwardMultiplicityBase.h"
#include "AliForwardCorrectionManager.h"
#include "AliForwardUtil.h"
#include "AliLog.h"
#include "AliAODHandler.h"
#include "AliInputEventHandler.h"
#include "AliAnalysisManager.h"
#include "AliFMDEventInspector.h"
#include "AliFMDSharingFilter.h"
#include "AliFMDDensityCalculator.h"
#include "AliFMDCorrector.h"
#include "AliFMDHistCollector.h"
#include "AliESDEvent.h"
#include <TROOT.h>
#include <TAxis.h>
#include <THStack.h>
#include <iostream>
#include <iomanip>

//====================================================================
AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name) 
  : AliAnalysisTaskSE(name), 
    fEnableLowFlux(false), 
    fFirstEvent(true),
    fCorrManager(0)
{
  // Set our persistent pointer 
  fCorrManager = &AliForwardCorrectionManager::Instance();
  fBranchNames = 
    "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
    "AliESDFMD.,SPDVertex.,PrimaryVertex.";
}

//____________________________________________________________________
AliForwardMultiplicityBase& 
AliForwardMultiplicityBase::operator=(const AliForwardMultiplicityBase& o)
{
  if (&o == this) return *this;
  fEnableLowFlux = o.fEnableLowFlux;
  fFirstEvent    = o.fFirstEvent;
  fCorrManager   = o.fCorrManager;
  return *this;
}

//____________________________________________________________________
Bool_t 
AliForwardMultiplicityBase::CheckCorrections(UInt_t what) const
{
  // 
  // Check if all needed corrections are there and accounted for.  If not,
  // do a Fatal exit 
  // 
  // Parameters:
  //    what Which corrections is needed
  // 
  // Return:
  //    true if all present, false otherwise
  //  

  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  // Check that we have the energy loss fits, needed by 
  //   AliFMDSharingFilter 
  //   AliFMDDensityCalculator 
  if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) { 
    AliFatal(Form("No energy loss fits"));
    return false;
  }
  // Check that we have the double hit correction - (optionally) used by 
  //  AliFMDDensityCalculator 
  if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit()) {
    AliFatal("No double hit corrections"); 
    return false;
  }
  // Check that we have the secondary maps, needed by 
  //   AliFMDCorrector 
  //   AliFMDHistCollector
  if (what & AliForwardCorrectionManager::kSecondaryMap && 
      !fcm.GetSecondaryMap()) {
    AliFatal("No secondary corrections");
    return false;
  }
  // Check that we have the vertex bias correction, needed by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kVertexBias && 
      !fcm.GetVertexBias()) { 
    AliFatal("No event vertex bias corrections");
    return false;
  }
  // Check that we have the merging efficiencies, optionally used by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kMergingEfficiency && 
      !fcm.GetMergingEfficiency()) {
    AliFatal("No merging efficiencies");
    return false;
  }
  // Check that we have the acceptance correction, needed by 
  //   AliFMDCorrector 
  if (what & AliForwardCorrectionManager::kAcceptance && 
      !fcm.GetAcceptance()) { 
    AliFatal("No acceptance corrections");
    return false;
  }
  return true;
}
//____________________________________________________________________
Bool_t
AliForwardMultiplicityBase::ReadCorrections(const TAxis*& pe, 
					    const TAxis*& pv, 
					    Bool_t        mc)
{
  //
  // Read corrections
  //
  //
  UInt_t what = AliForwardCorrectionManager::kAll;
  if (!fEnableLowFlux)
    what ^= AliForwardCorrectionManager::kDoubleHit;
  if (!GetCorrections().IsUseVertexBias())
    what ^= AliForwardCorrectionManager::kVertexBias;
  if (!GetCorrections().IsUseAcceptance())
    what ^= AliForwardCorrectionManager::kAcceptance;
  if (!GetCorrections().IsUseMergingEfficiency())
    what ^= AliForwardCorrectionManager::kMergingEfficiency;

  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
		GetEventInspector().GetEnergy(),
		GetEventInspector().GetField(),
		mc,
		what)) return false;
  if (!CheckCorrections(what)) return false;

  // Sett our persistency pointer 
  // fCorrManager = &fcm;

  // Get the eta axis from the secondary maps - if read in
  if (!pe) {
    pe = fcm.GetEtaAxis();
    if (!pe) AliFatal("No eta axis defined");
  }
  // Get the vertex axis from the secondary maps - if read in
  if (!pv) {
    pv = fcm.GetVertexAxis();
    if (!pv) AliFatal("No vertex axis defined");
  }

  return true;
}
//____________________________________________________________________
AliESDEvent*
AliForwardMultiplicityBase::GetESDEvent()
{
  //
  // 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;
  }

  // On the first event, initialize the parameters
  if (fFirstEvent && esd->GetESDRun()) {
    GetEventInspector().ReadRunDetails(esd);

    AliInfo(Form("Initializing with parameters from the ESD:\n"
                 "         AliESDEvent::GetBeamEnergy()   ->%f\n"
                 "         AliESDEvent::GetBeamType()     ->%s\n"
                 "         AliESDEvent::GetCurrentL3()    ->%f\n"
                 "         AliESDEvent::GetMagneticField()->%f\n"
                 "         AliESDEvent::GetRunNumber()    ->%d\n",
                 esd->GetBeamEnergy(),
                 esd->GetBeamType(),
                 esd->GetCurrentL3(),
                 esd->GetMagneticField(),
                 esd->GetRunNumber()));

    fFirstEvent = false;

    InitializeSubs();
  }
  return esd;
}
//____________________________________________________________________
void
AliForwardMultiplicityBase::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);
}

//____________________________________________________________________
void
AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input, 
					   const char*  inName,
					   TList*       output,
					   const char*  outName,
					   Int_t        style) const
{
  // Make dN/deta for each ring found in the input list.  
  // 
  // A stack of all the dN/deta is also made for easy drawing. 
  // 
  // Note, that the distributions are normalised to the number of
  // observed events only - they should be corrected for 
  if (!input) return;
  TList* list = static_cast<TList*>(input->FindObject(inName));
  if (!list) { 
    AliWarning(Form("No list %s found in %s", inName, input->GetName()));
    return;
  }
  
  TList* out = new TList;
  out->SetName(outName);
  out->SetOwner();
  output->Add(out);

  THStack*     dndetaRings = new THStack("all", "dN/d#eta per ring");
  const char*  names[]     = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
  const char** ptr         = names;
  
  while (*ptr) { 
    TList* thisList = new TList;
    thisList->SetOwner();
    thisList->SetName(*ptr);
    out->Add(thisList);

    TH2D* h = static_cast<TH2D*>(list->FindObject(Form("%s_cache", *ptr)));
    if (!h) { 
      AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
      ptr++;
      continue;
    }
    TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
    copy->SetDirectory(0);
    thisList->Add(copy);
    
    TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
    for (Int_t i = 1; i <= copy->GetNbinsX(); i++) { 
      for (Int_t j = 1; j <= copy->GetNbinsY(); j++) { 
	Double_t c = copy->GetBinContent(i, j);
	Double_t e = copy->GetBinError(i, j);
	Double_t a = norm->GetBinContent(i);
	copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
	copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
      }
    }

    TH1D* res  =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
						     h->GetNbinsY(),"e"));
    TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
    res->SetTitle(*ptr);
    res->Scale(1., "width");
    copy->Scale(1., "width");
    proj->Scale(1. / norm->GetMaximum(), "width");
    norm->Scale(1. / norm->GetMaximum());

    res->SetMarkerStyle(style);
    norm->SetDirectory(0);
    res->SetDirectory(0);
    proj->SetDirectory(0);
    thisList->Add(norm);
    thisList->Add(res);
    thisList->Add(proj);
    dndetaRings->Add(res);
    ptr++;
  }
  out->Add(dndetaRings);
}

//____________________________________________________________________
void
AliForwardMultiplicityBase::Print(Option_t* option) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
  
  std::cout << ClassName() << ": " << GetName() << "\n" 
	    << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no") 
	    << "\n"
	    << "  Off-line trigger mask:  0x" 
	    << std::hex     << std::setfill('0') 
	    << std::setw (8) << fOfflineTriggerMask 
	    << std::dec     << std::setfill (' ') << std::endl;
  gROOT->IncreaseDirLevel();
  if (fCorrManager) fCorrManager->Print();
  else  
    std::cout << "  Correction manager not set yet" << std::endl;
  GetEventInspector()   .Print(option);
  GetSharingFilter()    .Print(option);
  GetDensityCalculator().Print(option);
  GetCorrections()      .Print(option);
  GetHistCollector()    .Print(option);
  gROOT->DecreaseDirLevel();
}


//
// EOF
//
Commit	Line	Data
1a26066e	1	//====================================================================
7984e5f7	2	//
	3	// Base class for classes that calculate the multiplicity in the
	4	// forward regions event-by-event
	5	//
	6	// Inputs:
	7	// - AliESDEvent
	8	//
	9	// Outputs:
	10	// - AliAODForwardMult
	11	//
	12	// Histograms
	13	//
	14	// Corrections used
1a26066e	15	#include "AliForwardMultiplicityBase.h"
1174780f	16	#include "AliForwardCorrectionManager.h"
fb3430ac	17	#include "AliForwardUtil.h"
1a26066e	18	#include "AliLog.h"
	19	#include "AliAODHandler.h"
	20	#include "AliInputEventHandler.h"
	21	#include "AliAnalysisManager.h"
1174780f	22	#include "AliFMDEventInspector.h"
1174780f	23	#include "AliFMDSharingFilter.h"
1174780f	24	#include "AliFMDDensityCalculator.h"
72cc12cd	25	#include "AliFMDCorrector.h"
1174780f	26	#include "AliFMDHistCollector.h"
7ec4d843	27	#include "AliESDEvent.h"
1a26066e	28	#include <TROOT.h>
fb3430ac	29	#include <TAxis.h>
5bb5d1f6	30	#include <THStack.h>
3c5497d0	31	#include <iostream>
3c5497d0	32	#include <iomanip>
1a26066e	33
1a26066e	34	//====================================================================
19abe41d	35	AliForwardMultiplicityBase::AliForwardMultiplicityBase(const char* name)
19abe41d	36	: AliAnalysisTaskSE(name),
9d05ffeb	37	fEnableLowFlux(false),
19abe41d	38	fFirstEvent(true),
	39	fCorrManager(0)
	40	{
	41	// Set our persistent pointer
	42	fCorrManager = &AliForwardCorrectionManager::Instance();
331c3502	43	fBranchNames =
	44	"ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
	45	"AliESDFMD.,SPDVertex.,PrimaryVertex.";
19abe41d	46	}
19abe41d	47
fb3430ac	48	//____________________________________________________________________
	49	AliForwardMultiplicityBase&
	50	AliForwardMultiplicityBase::operator=(const AliForwardMultiplicityBase& o)
	51	{
d015ecfe	52	if (&o == this) return *this;
fb3430ac	53	fEnableLowFlux = o.fEnableLowFlux;
	54	fFirstEvent = o.fFirstEvent;
	55	fCorrManager = o.fCorrManager;
	56	return *this;
	57	}
	58
19abe41d	59	//____________________________________________________________________
1174780f	60	Bool_t
	61	AliForwardMultiplicityBase::CheckCorrections(UInt_t what) const
	62	{
	63	//
	64	// Check if all needed corrections are there and accounted for. If not,
	65	// do a Fatal exit
	66	//
	67	// Parameters:
	68	// what Which corrections is needed
	69	//
	70	// Return:
	71	// true if all present, false otherwise
	72	//
	73
	74	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	75	// Check that we have the energy loss fits, needed by
	76	// AliFMDSharingFilter
	77	// AliFMDDensityCalculator
	78	if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) {
	79	AliFatal(Form("No energy loss fits"));
	80	return false;
	81	}
	82	// Check that we have the double hit correction - (optionally) used by
	83	// AliFMDDensityCalculator
	84	if (what & AliForwardCorrectionManager::kDoubleHit && !fcm.GetDoubleHit()) {
	85	AliFatal("No double hit corrections");
	86	return false;
	87	}
	88	// Check that we have the secondary maps, needed by
72cc12cd	89	// AliFMDCorrector
1174780f	90	// AliFMDHistCollector
	91	if (what & AliForwardCorrectionManager::kSecondaryMap &&
	92	!fcm.GetSecondaryMap()) {
	93	AliFatal("No secondary corrections");
	94	return false;
	95	}
	96	// Check that we have the vertex bias correction, needed by
72cc12cd	97	// AliFMDCorrector
1174780f	98	if (what & AliForwardCorrectionManager::kVertexBias &&
	99	!fcm.GetVertexBias()) {
	100	AliFatal("No event vertex bias corrections");
	101	return false;
	102	}
	103	// Check that we have the merging efficiencies, optionally used by
72cc12cd	104	// AliFMDCorrector
1174780f	105	if (what & AliForwardCorrectionManager::kMergingEfficiency &&
	106	!fcm.GetMergingEfficiency()) {
	107	AliFatal("No merging efficiencies");
	108	return false;
	109	}
cc83fca2	110	// Check that we have the acceptance correction, needed by
	111	// AliFMDCorrector
	112	if (what & AliForwardCorrectionManager::kAcceptance &&
	113	!fcm.GetAcceptance()) {
	114	AliFatal("No acceptance corrections");
	115	return false;
	116	}
1174780f	117	return true;
1174780f	118	}
7ec4d843	119	//____________________________________________________________________
7ec4d843	120	Bool_t
19abe41d	121	AliForwardMultiplicityBase::ReadCorrections(const TAxis*& pe,
	122	const TAxis*& pv,
	123	Bool_t mc)
7ec4d843	124	{
fb3430ac	125	//
	126	// Read corrections
	127	//
	128	//
7ec4d843	129	UInt_t what = AliForwardCorrectionManager::kAll;
	130	if (!fEnableLowFlux)
	131	what ^= AliForwardCorrectionManager::kDoubleHit;
7ec4d843	132	if (!GetCorrections().IsUseVertexBias())
	133	what ^= AliForwardCorrectionManager::kVertexBias;
	134	if (!GetCorrections().IsUseAcceptance())
	135	what ^= AliForwardCorrectionManager::kAcceptance;
	136	if (!GetCorrections().IsUseMergingEfficiency())
	137	what ^= AliForwardCorrectionManager::kMergingEfficiency;
	138
	139	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	140	if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
19abe41d	141	GetEventInspector().GetEnergy(),
	142	GetEventInspector().GetField(),
	143	mc,
	144	what)) return false;
7ec4d843	145	if (!CheckCorrections(what)) return false;
7ec4d843	146
19abe41d	147	// Sett our persistency pointer
	148	// fCorrManager = &fcm;
	149
7ec4d843	150	// Get the eta axis from the secondary maps - if read in
	151	if (!pe) {
	152	pe = fcm.GetEtaAxis();
	153	if (!pe) AliFatal("No eta axis defined");
	154	}
	155	// Get the vertex axis from the secondary maps - if read in
	156	if (!pv) {
	157	pv = fcm.GetVertexAxis();
	158	if (!pv) AliFatal("No vertex axis defined");
	159	}
	160
	161	return true;
	162	}
	163	//____________________________________________________________________
	164	AliESDEvent*
	165	AliForwardMultiplicityBase::GetESDEvent()
	166	{
fb3430ac	167	//
	168	// Get the ESD event. IF this is the first event, initialise
	169	//
7ec4d843	170	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	171	if (!esd) {
	172	AliWarning("No ESD event found for input event");
	173	return 0;
	174	}
1174780f	175
7ec4d843	176	// On the first event, initialize the parameters
	177	if (fFirstEvent && esd->GetESDRun()) {
	178	GetEventInspector().ReadRunDetails(esd);
	179
	180	AliInfo(Form("Initializing with parameters from the ESD:\n"
	181	" AliESDEvent::GetBeamEnergy() ->%f\n"
	182	" AliESDEvent::GetBeamType() ->%s\n"
	183	" AliESDEvent::GetCurrentL3() ->%f\n"
	184	" AliESDEvent::GetMagneticField()->%f\n"
	185	" AliESDEvent::GetRunNumber() ->%d\n",
	186	esd->GetBeamEnergy(),
	187	esd->GetBeamType(),
	188	esd->GetCurrentL3(),
	189	esd->GetMagneticField(),
	190	esd->GetRunNumber()));
	191
	192	fFirstEvent = false;
	193
	194	InitializeSubs();
	195	}
	196	return esd;
	197	}
1174780f	198	//____________________________________________________________________
1a26066e	199	void
	200	AliForwardMultiplicityBase::MarkEventForStore() const
	201	{
	202	// Make sure the AOD tree is filled
	203	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	204	AliAODHandler* ah =
	205	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	206	if (!ah)
	207	AliFatal("No AOD output handler set in analysis manager");
	208
	209	ah->SetFillAOD(kTRUE);
	210	}
	211
5bb5d1f6	212	//____________________________________________________________________
	213	void
	214	AliForwardMultiplicityBase::MakeRingdNdeta(const TList* input,
	215	const char* inName,
	216	TList* output,
	217	const char* outName,
	218	Int_t style) const
	219	{
	220	// Make dN/deta for each ring found in the input list.
	221	//
	222	// A stack of all the dN/deta is also made for easy drawing.
	223	//
	224	// Note, that the distributions are normalised to the number of
	225	// observed events only - they should be corrected for
	226	if (!input) return;
	227	TList* list = static_cast<TList*>(input->FindObject(inName));
	228	if (!list) {
	229	AliWarning(Form("No list %s found in %s", inName, input->GetName()));
	230	return;
	231	}
	232
	233	TList* out = new TList;
	234	out->SetName(outName);
	235	out->SetOwner();
	236	output->Add(out);
	237
	238	THStack* dndetaRings = new THStack("all", "dN/d#eta per ring");
	239	const char* names[] = { "FMD1I", "FMD2I", "FMD2O", "FMD3I", "FMD3O", 0 };
	240	const char** ptr = names;
	241
	242	while (*ptr) {
	243	TList* thisList = new TList;
	244	thisList->SetOwner();
	245	thisList->SetName(*ptr);
	246	out->Add(thisList);
	247
	248	TH2D* h = static_cast<TH2D>(list->FindObject(Form("%s_cache", ptr)));
	249	if (!h) {
	250	AliWarning(Form("Didn't find %s_cache in %s", *ptr, list->GetName()));
	251	ptr++;
	252	continue;
	253	}
	254	TH2D* copy = static_cast<TH2D*>(h->Clone("sum"));
	255	copy->SetDirectory(0);
	256	thisList->Add(copy);
	257
	258	TH1D* norm =static_cast<TH1D*>(h->ProjectionX("norm", 0, 0, ""));
	259	for (Int_t i = 1; i <= copy->GetNbinsX(); i++) {
	260	for (Int_t j = 1; j <= copy->GetNbinsY(); j++) {
	261	Double_t c = copy->GetBinContent(i, j);
	262	Double_t e = copy->GetBinError(i, j);
	263	Double_t a = norm->GetBinContent(i);
	264	copy->SetBinContent(i, j, a <= 0 ? 0 : c / a);
	265	copy->SetBinError(i, j, a <= 0 ? 0 : e / a);
	266	}
	267	}
	268
	269	TH1D* res =static_cast<TH1D*>(copy->ProjectionX("dndeta",1,
	270	h->GetNbinsY(),"e"));
	271	TH1D* proj =static_cast<TH1D*>(h->ProjectionX("proj",1,h->GetNbinsY(),"e"));
	272	res->SetTitle(*ptr);
	273	res->Scale(1., "width");
	274	copy->Scale(1., "width");
	275	proj->Scale(1. / norm->GetMaximum(), "width");
276	norm->Scale(1. / norm->GetMaximum());
277
278	res->SetMarkerStyle(style);
279	norm->SetDirectory(0);
280	res->SetDirectory(0);
281	proj->SetDirectory(0);
282	thisList->Add(norm);
283	thisList->Add(res);
284	thisList->Add(proj);
285	dndetaRings->Add(res);
286	ptr++;
287	}
288	out->Add(dndetaRings);
289	}
290
1a26066e	291	//____________________________________________________________________
	292	void
	293	AliForwardMultiplicityBase::Print(Option_t* option) const
	294	{
7984e5f7	295	//
	296	// Print information
	297	//
	298	// Parameters:
	299	// option Not used
	300	//
3c5497d0	301
1f480471	302	std::cout << ClassName() << ": " << GetName() << "\n"
9bff843b	303	<< " Enable low flux code: " << (fEnableLowFlux ? "yes" : "no")
9bff843b	304	<< "\n"
3c5497d0	305	<< " Off-line trigger mask: 0x"
	306	<< std::hex << std::setfill('0')
	307	<< std::setw (8) << fOfflineTriggerMask
	308	<< std::dec << std::setfill (' ') << std::endl;
1174780f	309	gROOT->IncreaseDirLevel();
3c5497d0	310	if (fCorrManager) fCorrManager->Print();
9bff843b	311	else
9bff843b	312	std::cout << " Correction manager not set yet" << std::endl;
1174780f	313	GetEventInspector() .Print(option);
1174780f	314	GetSharingFilter() .Print(option);
	315	GetDensityCalculator().Print(option);
	316	GetCorrections() .Print(option);
	317	GetHistCollector() .Print(option);
	318	gROOT->DecreaseDirLevel();
1a26066e	319	}
1a26066e	320
1174780f	321
1a26066e	322	//
	323	// EOF
	324	//