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

// 
// Histogram and fit the energy loss distributions for the FMD
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - None
// 
// Histograms:
//   
// Corrections used:
//   - None
// 
// 
//
#include "AliFMDEnergyFitterTask.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliAODForwardMult.h"
#include "AliForwardCorrectionManager.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisDataContainer.h"
#include <TH1.h>
#include <TDirectory.h>
#include <TTree.h>
#include <TFile.h>
#include "AliMCEvent.h"
#include "AliGenHijingEventHeader.h"
#include "AliHeader.h"
#include <iostream>

//====================================================================
AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
  : AliAnalysisTaskSE(),
    fFirstEvent(true),
    fEventInspector(),
    fEnergyFitter(),
    fList(0),
    fbLow(0),
    fbHigh(100)
{
  // 
  // Constructor
  //
  DGUARD(fDebug,0,"Default CTOR of AliFMDEnergyFitterTask");
}

//____________________________________________________________________
AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name)
  : AliAnalysisTaskSE(name), 
    fFirstEvent(true),
    fEventInspector("event"),
    fEnergyFitter("energy"),
    fList(0),
    fbLow(0),
    fbHigh(100)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug,0,"Named CTOR of AliFMDEnergyFitterTask: %s", name);
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  fBranchNames = 
    "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
    "AliESDFMD.,SPDVertex.,PrimaryVertex.";
}

//____________________________________________________________________
AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o)
  : AliAnalysisTaskSE(o),
    fFirstEvent(o.fFirstEvent),
    fEventInspector(o.fEventInspector),
    fEnergyFitter(o.fEnergyFitter),
    fList(o.fList),
    fbLow(o.fbLow),
    fbHigh(o.fbHigh)
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
  DGUARD(fDebug,0,"COPY CTOR of AliFMDEnergyFitterTask");
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
}

//____________________________________________________________________
AliFMDEnergyFitterTask&
AliFMDEnergyFitterTask::operator=(const AliFMDEnergyFitterTask& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object 
  //
  DGUARD(fDebug,3,"Assignment of AliFMDEnergyFitterTask");
  if (&o == this) return *this; 
  AliAnalysisTaskSE::operator=(o);

  fFirstEvent        = o.fFirstEvent;
  fEventInspector    = o.fEventInspector;
  fEnergyFitter      = o.fEnergyFitter;
  fList              = o.fList;
  fbLow              = o.fbLow;
  fbHigh             = o.fbHigh;

  return *this;
}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
{
  // 
  // Set the debug level 
  // 
  // Parameters:
  //    dbg Debug level
  //
  fEventInspector.SetDebug(dbg);
  fEnergyFitter.SetDebug(dbg);
}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::Init()
{
  // 
  // Initialize the task 
  // 
  //
  DGUARD(fDebug,1,"Initialize of AliFMDEnergyFitterTask");
  fFirstEvent = true;
}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::InitializeSubs()
{
  // 
  // Initialise the sub objects and stuff.  Called on first event 
  // 
  //
  DGUARD(fDebug,1,"Initialize subs of AliFMDEnergyFitterTask");
  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  UShort_t sys = fEventInspector.GetCollisionSystem();
  UShort_t sNN = fEventInspector.GetEnergy();
  Short_t  fld = fEventInspector.GetField();
  fcm.Init(sys, sNN, fld, 0);
  TAxis eAxis(0,0,0);
  TAxis vAxis(10,-10,10);
  fEnergyFitter.Init(eAxis);
  fEventInspector.Init(vAxis);

}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::UserCreateOutputObjects()
{
  // 
  // Create output objects 
  // 
  //
  DGUARD(fDebug,1,"Create output objects of AliFMDEnergyFitterTask");
  fList = new TList;
  fList->SetOwner();

  fEventInspector.DefineOutput(fList);
  fEnergyFitter.DefineOutput(fList);

  PostData(1, fList);
}
//____________________________________________________________________
void
AliFMDEnergyFitterTask::UserExec(Option_t*)
{
  // 
  // Process each event 
  // 
  // Parameters:
  //    option Not used
  //  

  // static Int_t cnt = 0;
  // cnt++;
  // Get the input data 
  DGUARD(fDebug,1,"Analyse event of AliFMDEnergyFitterTask");
  
  AliMCEvent* mcevent = MCEvent();
  if(mcevent) {
    AliHeader* header            = mcevent->Header();
    AliGenHijingEventHeader* hijingHeader = 
      dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
    if(hijingHeader) {
      Float_t b = hijingHeader->ImpactParameter();
      if(b<fbLow || b>fbHigh) return;
      else
	std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
    }
    
  }
  
  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
  // AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
  if (!esd) { 
    AliWarning("No ESD event found for input event");
    return;
  }

  // On the first event, initialize the parameters 
  if (fFirstEvent && esd->GetESDRun()) { 
    fEventInspector.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()));

	      

    // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
    // pars->SetParametersFromESD(esd);
    // pars->PrintStatus();
    fFirstEvent = false;

    InitializeSubs();
  }
  Bool_t   lowFlux   = kFALSE;
  UInt_t   triggers  = 0;
  UShort_t ivz       = 0;
  Double_t vz        = 0;
  Double_t cent      = 0;
  UShort_t nClusters = 0;
  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
					       ivz, vz, cent, nClusters);
  if (found & AliFMDEventInspector::kNoEvent)    return;
  if (found & AliFMDEventInspector::kNoTriggers) return;
  if (found & AliFMDEventInspector::kNoSPD)     return;
  if (found & AliFMDEventInspector::kNoFMD)     return;
  if (found & AliFMDEventInspector::kNoVertex)  return;
  if (found & AliFMDEventInspector::kBadVertex) return;
  
  //  if(cent > 0) {
  //  if( cent < 40 || cent >50 ) return;
  //  else std::cout<<"selecting event with cent "<<cent<<std::endl;
  // }
  
  // Get FMD data 
  AliESDFMD* esdFMD = esd->GetFMDData();
  // Do the energy stuff 
  if (!fEnergyFitter.Accumulate(*esdFMD, cent, 
				triggers & AliAODForwardMult::kEmpty)){
    AliWarning("Energy fitter failed");
    return;
  }
  PostData(1, fList);
}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::Terminate(Option_t*)
{
  // 
  // End of job
  // 
  // Parameters:
  //    option Not used 
  //
  DGUARD(fDebug,1,"Processing merged output of AliFMDEnergyFitterTask");

  TList* list = dynamic_cast<TList*>(GetOutputData(1));
  if (!list) {
    AliError(Form("No output list defined (%p)", GetOutputData(1)));
    if (GetOutputData(1)) GetOutputData(1)->Print();
    return;
  }
  
  AliInfo("Fitting energy loss spectra");
  fEnergyFitter.Fit(list);

  // Make a deep copy and post that as output 2 
  TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults", 
						      list->GetName())));
  list2->SetOwner();
  PostData(2, list2);
}

//____________________________________________________________________
void
AliFMDEnergyFitterTask::Print(Option_t*) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
}

//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// Histogram and fit the energy loss distributions for the FMD
	3	//
	4	// Inputs:
	5	// - AliESDEvent
	6	//
	7	// Outputs:
	8	// - None
	9	//
	10	// Histograms:
	11	//
	12	// Corrections used:
	13	// - None
	14	//
	15	//
	16	//
0bd4b00f	17	#include "AliFMDEnergyFitterTask.h"
0bd4b00f	18	#include "AliLog.h"
0bd4b00f	19	#include "AliESDEvent.h"
	20	#include "AliAODForwardMult.h"
	21	#include "AliForwardCorrectionManager.h"
	22	#include "AliAnalysisManager.h"
72cc12cd	23	#include "AliAnalysisDataSlot.h"
72cc12cd	24	#include "AliAnalysisDataContainer.h"
0bd4b00f	25	#include <TH1.h>
	26	#include <TDirectory.h>
	27	#include <TTree.h>
72cc12cd	28	#include <TFile.h>
e58000b7	29	#include "AliMCEvent.h"
	30	#include "AliGenHijingEventHeader.h"
	31	#include "AliHeader.h"
	32	#include <iostream>
0bd4b00f	33
	34	//====================================================================
	35	AliFMDEnergyFitterTask::AliFMDEnergyFitterTask()
	36	: AliAnalysisTaskSE(),
	37	fFirstEvent(true),
	38	fEventInspector(),
	39	fEnergyFitter(),
6f4c294b	40	fList(0),
	41	fbLow(0),
	42	fbHigh(100)
0bd4b00f	43	{
7984e5f7	44	//
	45	// Constructor
	46	//
6ab100ec	47	DGUARD(fDebug,0,"Default CTOR of AliFMDEnergyFitterTask");
0bd4b00f	48	}
	49
	50	//____________________________________________________________________
	51	AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const char* name)
	52	: AliAnalysisTaskSE(name),
	53	fFirstEvent(true),
	54	fEventInspector("event"),
	55	fEnergyFitter("energy"),
e58000b7	56	fList(0),
	57	fbLow(0),
	58	fbHigh(100)
0bd4b00f	59	{
7984e5f7	60	//
	61	// Constructor
	62	//
	63	// Parameters:
	64	// name Name of task
	65	//
6ab100ec	66	DGUARD(fDebug,0,"Named CTOR of AliFMDEnergyFitterTask: %s", name);
0bd4b00f	67	DefineOutput(1, TList::Class());
8e2bb72a	68	DefineOutput(2, TList::Class());
5db7b981	69	fBranchNames =
	70	"ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
	71	"AliESDFMD.,SPDVertex.,PrimaryVertex.";
0bd4b00f	72	}
	73
	74	//____________________________________________________________________
	75	AliFMDEnergyFitterTask::AliFMDEnergyFitterTask(const AliFMDEnergyFitterTask& o)
	76	: AliAnalysisTaskSE(o),
	77	fFirstEvent(o.fFirstEvent),
	78	fEventInspector(o.fEventInspector),
	79	fEnergyFitter(o.fEnergyFitter),
6f4c294b	80	fList(o.fList),
	81	fbLow(o.fbLow),
	82	fbHigh(o.fbHigh)
0bd4b00f	83	{
7984e5f7	84	//
	85	// Copy constructor
	86	//
	87	// Parameters:
	88	// o Object to copy from
	89	//
6ab100ec	90	DGUARD(fDebug,0,"COPY CTOR of AliFMDEnergyFitterTask");
0bd4b00f	91	DefineOutput(1, TList::Class());
8e2bb72a	92	DefineOutput(2, TList::Class());
0bd4b00f	93	}
	94
	95	//____________________________________________________________________
	96	AliFMDEnergyFitterTask&
	97	AliFMDEnergyFitterTask::operator=(const AliFMDEnergyFitterTask& o)
	98	{
7984e5f7	99	//
	100	// Assignment operator
	101	//
	102	// Parameters:
	103	// o Object to assign from
	104	//
	105	// Return:
	106	// Reference to this object
	107	//
6ab100ec	108	DGUARD(fDebug,3,"Assignment of AliFMDEnergyFitterTask");
d015ecfe	109	if (&o == this) return *this;
0bd4b00f	110	AliAnalysisTaskSE::operator=(o);
	111
	112	fFirstEvent = o.fFirstEvent;
	113	fEventInspector = o.fEventInspector;
	114	fEnergyFitter = o.fEnergyFitter;
	115	fList = o.fList;
6f4c294b	116	fbLow = o.fbLow;
6f4c294b	117	fbHigh = o.fbHigh;
0bd4b00f	118
	119	return *this;
	120	}
	121
	122	//____________________________________________________________________
	123	void
	124	AliFMDEnergyFitterTask::SetDebug(Int_t dbg)
	125	{
7984e5f7	126	//
	127	// Set the debug level
	128	//
	129	// Parameters:
	130	// dbg Debug level
	131	//
0bd4b00f	132	fEventInspector.SetDebug(dbg);
	133	fEnergyFitter.SetDebug(dbg);
	134	}
	135
	136	//____________________________________________________________________
	137	void
	138	AliFMDEnergyFitterTask::Init()
	139	{
7984e5f7	140	//
	141	// Initialize the task
	142	//
	143	//
6ab100ec	144	DGUARD(fDebug,1,"Initialize of AliFMDEnergyFitterTask");
0bd4b00f	145	fFirstEvent = true;
	146	}
	147
	148	//____________________________________________________________________
	149	void
	150	AliFMDEnergyFitterTask::InitializeSubs()
	151	{
7984e5f7	152	//
	153	// Initialise the sub objects and stuff. Called on first event
	154	//
	155	//
6ab100ec	156	DGUARD(fDebug,1,"Initialize subs of AliFMDEnergyFitterTask");
0bd4b00f	157	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
f7cfc454	158	UShort_t sys = fEventInspector.GetCollisionSystem();
	159	UShort_t sNN = fEventInspector.GetEnergy();
	160	Short_t fld = fEventInspector.GetField();
	161	fcm.Init(sys, sNN, fld, 0);
0bd4b00f	162	TAxis eAxis(0,0,0);
	163	TAxis vAxis(10,-10,10);
	164	fEnergyFitter.Init(eAxis);
	165	fEventInspector.Init(vAxis);
f7cfc454	166
0bd4b00f	167	}
	168
	169	//____________________________________________________________________
	170	void
	171	AliFMDEnergyFitterTask::UserCreateOutputObjects()
	172	{
7984e5f7	173	//
	174	// Create output objects
	175	//
	176	//
6ab100ec	177	DGUARD(fDebug,1,"Create output objects of AliFMDEnergyFitterTask");
0bd4b00f	178	fList = new TList;
9b07e44d	179	fList->SetOwner();
0bd4b00f	180
	181	fEventInspector.DefineOutput(fList);
	182	fEnergyFitter.DefineOutput(fList);
	183
	184	PostData(1, fList);
	185	}
	186	//____________________________________________________________________
	187	void
	188	AliFMDEnergyFitterTask::UserExec(Option_t*)
	189	{
7984e5f7	190	//
	191	// Process each event
	192	//
	193	// Parameters:
	194	// option Not used
	195	//
	196
0bd4b00f	197	// static Int_t cnt = 0;
	198	// cnt++;
	199	// Get the input data
6ab100ec	200	DGUARD(fDebug,1,"Analyse event of AliFMDEnergyFitterTask");
e58000b7	201
	202	AliMCEvent* mcevent = MCEvent();
	203	if(mcevent) {
	204	AliHeader* header = mcevent->Header();
5bb5d1f6	205	AliGenHijingEventHeader* hijingHeader =
5bb5d1f6	206	dynamic_cast<AliGenHijingEventHeader*>(header->GenEventHeader());
e58000b7	207	if(hijingHeader) {
	208	Float_t b = hijingHeader->ImpactParameter();
	209	if(b<fbLow \|\| b>fbHigh) return;
	210	else
	211	std::cout<<"Selecting event with impact parameter "<<b<<std::endl;
	212	}
	213
	214	}
	215
0bd4b00f	216	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	217	// AliInfo(Form("Event # %6d (esd=%p)", cnt, esd));
	218	if (!esd) {
	219	AliWarning("No ESD event found for input event");
	220	return;
	221	}
	222
	223	// On the first event, initialize the parameters
	224	if (fFirstEvent && esd->GetESDRun()) {
	225	fEventInspector.ReadRunDetails(esd);
	226
	227	AliInfo(Form("Initializing with parameters from the ESD:\n"
	228	" AliESDEvent::GetBeamEnergy() ->%f\n"
	229	" AliESDEvent::GetBeamType() ->%s\n"
	230	" AliESDEvent::GetCurrentL3() ->%f\n"
	231	" AliESDEvent::GetMagneticField()->%f\n"
	232	" AliESDEvent::GetRunNumber() ->%d\n",
	233	esd->GetBeamEnergy(),
	234	esd->GetBeamType(),
	235	esd->GetCurrentL3(),
	236	esd->GetMagneticField(),
	237	esd->GetRunNumber()));
	238
	239
	240
	241	// AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	242	// pars->SetParametersFromESD(esd);
	243	// pars->PrintStatus();
	244	fFirstEvent = false;
	245
	246	InitializeSubs();
	247	}
5bb5d1f6	248	Bool_t lowFlux = kFALSE;
	249	UInt_t triggers = 0;
	250	UShort_t ivz = 0;
	251	Double_t vz = 0;
	252	Double_t cent = 0;
	253	UShort_t nClusters = 0;
	254	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
	255	ivz, vz, cent, nClusters);
0bd4b00f	256	if (found & AliFMDEventInspector::kNoEvent) return;
	257	if (found & AliFMDEventInspector::kNoTriggers) return;
	258	if (found & AliFMDEventInspector::kNoSPD) return;
	259	if (found & AliFMDEventInspector::kNoFMD) return;
	260	if (found & AliFMDEventInspector::kNoVertex) return;
	261	if (found & AliFMDEventInspector::kBadVertex) return;
12fffad7	262
	263	// if(cent > 0) {
	264	// if( cent < 40 \|\| cent >50 ) return;
	265	// else std::cout<<"selecting event with cent "<<cent<<std::endl;
	266	// }
	267
0bd4b00f	268	// Get FMD data
	269	AliESDFMD* esdFMD = esd->GetFMDData();
	270	// Do the energy stuff
5e4d905e	271	if (!fEnergyFitter.Accumulate(*esdFMD, cent,
5e4d905e	272	triggers & AliAODForwardMult::kEmpty)){
0bd4b00f	273	AliWarning("Energy fitter failed");
	274	return;
	275	}
	276	PostData(1, fList);
	277	}
	278
	279	//____________________________________________________________________
	280	void
	281	AliFMDEnergyFitterTask::Terminate(Option_t*)
	282	{
7984e5f7	283	//
	284	// End of job
	285	//
	286	// Parameters:
	287	// option Not used
	288	//
6ab100ec	289	DGUARD(fDebug,1,"Processing merged output of AliFMDEnergyFitterTask");
6ab100ec	290
0bd4b00f	291	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	292	if (!list) {
	293	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	294	if (GetOutputData(1)) GetOutputData(1)->Print();
	295	return;
	296	}
	297
1c762251	298	AliInfo("Fitting energy loss spectra");
0bd4b00f	299	fEnergyFitter.Fit(list);
72cc12cd	300
8e2bb72a	301	// Make a deep copy and post that as output 2
	302	TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
	303	list->GetName())));
9b07e44d	304	list2->SetOwner();
8e2bb72a	305	PostData(2, list2);
0bd4b00f	306	}
	307
	308	//____________________________________________________________________
	309	void
	310	AliFMDEnergyFitterTask::Print(Option_t*) const
	311	{
7984e5f7	312	//
	313	// Print information
	314	//
	315	// Parameters:
	316	// option Not used
	317	//
0bd4b00f	318	}
	319
	320	//
	321	// EOF
	322	//