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

// 
// Calculate the multiplicity in the forward regions event-by-event 
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODForwardMult 
// 
// Histograms 
//   
// Corrections used 
//
#include "AliForwardQATask.h"
#include "AliForwardUtil.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
#include "AliForwardCorrectionManager.h"
#include "AliAnalysisManager.h"
#include "AliAODForwardMult.h"
#include <TH1.h>
#include <TDirectory.h>
#include <TTree.h>
#include <TROOT.h>
#include <TStopwatch.h>

//====================================================================
AliForwardQATask::AliForwardQATask()
  : AliAnalysisTaskSE(),
    fEnableLowFlux(false), 
    fFirstEvent(true),
    fCorrManager(0),
    fESDFMD(),
    fHistos(),
    fEventInspector(),
    fEnergyFitter(),
    fSharingFilter(),
    fDensityCalculator(),
    fList(0),
    fDebug(0)
{
  // 
  // Constructor
  //
}

//____________________________________________________________________
AliForwardQATask::AliForwardQATask(const char* name)
  : AliAnalysisTaskSE(name),
    fEnableLowFlux(false), 
    fFirstEvent(true),
    fCorrManager(0),
    fESDFMD(),
    fHistos(),
    fEventInspector("event"),
    fEnergyFitter("energy"),
    fSharingFilter("sharing"), 
    fDensityCalculator("density"),
    fList(0),
    fDebug(0)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  fCorrManager = &AliForwardCorrectionManager::Instance(); 
  fEnergyFitter.SetNParticles(1); // Just find the 1st peak 
  fEnergyFitter.SetDoMakeObject(false); 
  fEnergyFitter.SetUseIncreasingBins(true);
  fEnergyFitter.SetDoFits(kTRUE);
  fEnergyFitter.SetLowCut(0.4);
  fEnergyFitter.SetFitRangeBinWidth(4);
  fEnergyFitter.SetMinEntries(1000);
}

//____________________________________________________________________
AliForwardQATask::AliForwardQATask(const AliForwardQATask& o)
  : AliAnalysisTaskSE(o),
    fEnableLowFlux(o.fEnableLowFlux), 
    fFirstEvent(o.fFirstEvent),
    fCorrManager(o.fCorrManager),
    fESDFMD(o.fESDFMD),
    fHistos(o.fHistos),
    fEventInspector(o.fEventInspector),
    fEnergyFitter(o.fEnergyFitter),
    fSharingFilter(o.fSharingFilter),
    fDensityCalculator(o.fDensityCalculator),
    fList(o.fList),
    fDebug(o.fDebug) 
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
}

//____________________________________________________________________
AliForwardQATask&
AliForwardQATask::operator=(const AliForwardQATask& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object 
  //
  if (&o == this) return *this;
  AliAnalysisTaskSE::operator=(o);

  fEnableLowFlux     = o.fEnableLowFlux;
  fFirstEvent        = o.fFirstEvent;
  fCorrManager       = o.fCorrManager;
  fEventInspector    = o.fEventInspector;
  fEnergyFitter      = o.fEnergyFitter;
  fSharingFilter     = o.fSharingFilter;
  fDensityCalculator = o.fDensityCalculator;
  fHistos            = o.fHistos;
  fList              = o.fList;
  fDebug             = o.fDebug;

  return *this;
}

//____________________________________________________________________
void
AliForwardQATask::SetDebug(Int_t dbg)
{
  // 
  // Set debug level 
  // 
  // Parameters:
  //    dbg Debug level
  //
  fDebug = dbg;
  fEventInspector.SetDebug(dbg);
  fEnergyFitter.SetDebug(dbg);
  fSharingFilter.SetDebug(dbg);
  fDensityCalculator.SetDebug(dbg);
}

//____________________________________________________________________
Bool_t 
AliForwardQATask::CheckCorrections(UInt_t what)
{
  // 
  // 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) {
    if (!fcm.GetELossFit()) { 
      AliWarning("No energy loss fits");
      
      // Fall-back values if we do not have the energy loss fits 
      AliFMDMultCuts& sfLCuts = GetSharingFilter().GetLCuts();
      if (sfLCuts.GetMethod() != AliFMDMultCuts::kFixed) { 
	Double_t cut = 0.3;
	AliWarningF("Using fixed cut @ %f for the lower bound "
		    "of the sharing filter", cut);
	sfLCuts.SetMultCuts(cut);
      }
      AliFMDMultCuts& sfHCuts = GetSharingFilter().GetHCuts();
      if (sfHCuts.GetMethod() != AliFMDMultCuts::kFixed) { 
	Double_t cut = 100;
	AliWarningF("Using fixed cut @ %f for the upper bound "
		    "of the sharing filter", cut);
	sfHCuts.SetMultCuts(cut);
      }
      AliFMDMultCuts& dcCuts  = GetDensityCalculator().GetCuts();
      if (dcCuts.GetMethod() != AliFMDMultCuts::kFixed) { 
	Double_t cut = 0.3;
	AliWarningF("Using fixed cut @ %f for the lower bound "
		    "of the density calculator", cut);
	dcCuts.SetMultCuts(cut);
      }
    }
    else 
      fcm.GetELossFit()->CacheBins(GetDensityCalculator().GetMinQuality());
  }

  return true;
}

//____________________________________________________________________
Bool_t
AliForwardQATask::ReadCorrections(const TAxis*& pe, 
				  const TAxis*& pv, 
				  Bool_t        mc,
				  Bool_t        sat)
{
  //
  // Read corrections
  //
  //
  UInt_t what = AliForwardCorrectionManager::kAll;
  what ^= AliForwardCorrectionManager::kDoubleHit;
  what ^= AliForwardCorrectionManager::kVertexBias;
  what ^= AliForwardCorrectionManager::kAcceptance;
  what ^= AliForwardCorrectionManager::kMergingEfficiency;

  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  // We allow fall-back queries so that we may proceed in case we have no 
  // valid corrections 
  fcm.SetEnableFallBack(true); 
  if (!fcm.Init(GetEventInspector().GetRunNumber(), 
		GetEventInspector().GetCollisionSystem(),
		GetEventInspector().GetEnergy(),
		GetEventInspector().GetField(),
		mc,
		sat,
		what,
		false)) 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*
AliForwardQATask::GetESDEvent()
{
  //
  // Get the ESD event. IF this is the first event, initialise
  //
  DGUARD(fDebug,2,"Get the ESD event");
  if (IsZombie()) {
    DMSG(fDebug,3,"We're a Zombie - bailing out");
    return 0;
  }
  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);

    AliInfoF("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",
	     esd->GetBeamEnergy(),
	     esd->GetBeamType(),
	     esd->GetCurrentL3(),
	     esd->GetMagneticField(),
	     esd->GetRunNumber());


    if (!SetupForData()) {
      AliWarning("Initialisation of sub algorithms failed!");
      SetZombie(true);
      esd = 0;
      return 0;
    }
    AliInfoF("Clearing first event flag from %s to false", 
	     fFirstEvent ? "true" : "false");
    fFirstEvent = false;
  }
  return esd;
}
//____________________________________________________________________
Bool_t
AliForwardQATask::SetupForData()
{
  // 
  // Initialise the sub objects and stuff.  Called on first event 
  // 
  //
  const TAxis* pe = 0;
  const TAxis* pv = 0;

  if (!ReadCorrections(pe,pv))  { 
    AliWarning("Using default binning");
    pv = new TAxis(10,-10, 10);
    pe = new TAxis(240,-6,6);
  }

  fHistos.Init(*pe);

  fEventInspector.SetupForData(*pv);
  fEnergyFitter.SetupForData(*pe);
  fSharingFilter.SetupForData(*pe);
  fDensityCalculator.SetupForData(*pe);

  this->Print();

  return true;
}

//____________________________________________________________________
void
AliForwardQATask::UserCreateOutputObjects()
{
  // 
  // Create output objects 
  // 
  //
  fList = new TList;
  fList->SetOwner();
  
  fEventInspector.CreateOutputObjects(fList);
  fEnergyFitter.CreateOutputObjects(fList);
  fSharingFilter.CreateOutputObjects(fList);
  fDensityCalculator.CreateOutputObjects(fList);

  PostData(1, fList);
}
//____________________________________________________________________
void
AliForwardQATask::UserExec(Option_t*)
{
  // 
  // Process each event 
  // 
  // Parameters:
  //    option Not used
  //  
  DGUARD(fDebug,1,"Process the input event");

  // static Int_t cnt = 0;
  // cnt++;
  // Get the input data 
  AliESDEvent* esd = GetESDEvent();
  if (!esd) { 
    AliWarning("Got no ESD event");
    return;
  }
  if (fFirstEvent) { 
    // If the first event flag wasn't cleared in the above call to
    // GetESDEvent, we should not do anything, since nothing has been
    // initialised yet, so we opt out here (with a warning) 
    AliWarning("Nothing has been initialized yet, opt'ing out");
    return;
  }

  // Clear stuff 
  fHistos.Clear();
  fESDFMD.Clear();
  
  Bool_t   lowFlux   = kFALSE;
  UInt_t   triggers  = 0;
  UShort_t ivz       = 0;
  TVector3 ip;
  Double_t cent      = -1;
  UShort_t nClusters = 0;
  UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
					       ivz, ip, cent, nClusters);
  
  Bool_t ok = true;
  if (found & AliFMDEventInspector::kNoEvent)    ok = false;
  if (found & AliFMDEventInspector::kNoTriggers) ok = false;
  if (found & AliFMDEventInspector::kNoSPD)      ok = false;
  if (found & AliFMDEventInspector::kNoFMD)      ok = false;
  if (found & AliFMDEventInspector::kNoVertex)   ok = false;
  if (triggers & AliAODForwardMult::kPileUp)     ok = false;
  if (triggers & AliAODForwardMult::kA)          ok = false;
  if (triggers & AliAODForwardMult::kC)          ok = false;
  if (triggers & AliAODForwardMult::kE)          ok = false;
  if (!(triggers & AliAODForwardMult::kOffline)) ok = false;
  if (found & AliFMDEventInspector::kBadVertex)  ok = false;
  if (!ok) { 
    DMSG(fDebug,2,"Event failed selection: %s", 
	 AliFMDEventInspector::CodeString(found));
    return;
  }
  DMSG(fDebug,2,"Event triggers: %s", 
       AliAODForwardMult::GetTriggerString(triggers));

  // We we do not want to use low flux specific code, we disable it here. 
  if (!fEnableLowFlux) lowFlux = false;

  // Get FMD data 
  AliESDFMD* esdFMD = esd->GetFMDData();
  
  // Run the energy loss fitter 
  if (!fEnergyFitter.Accumulate(*esdFMD, cent, 
				triggers & AliAODForwardMult::kEmpty)) {
    AliWarning("Energy fitter failed");
    return;
  }
  
  //  // Apply the sharing filter (or hit merging or clustering if you like)
  if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) { 
    AliWarning("Sharing filter failed!");
    return;
  }
 
  // Calculate the inclusive charged particle density 
  if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) { 
    // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { 
    AliWarning("Density calculator failed!");
    return;
  }
  
  PostData(1, fList);
}

//____________________________________________________________________
void
AliForwardQATask::Terminate(Option_t*)
{
  // 
  // End of job
  // 
  // Parameters:
  //    option Not used 
  //
  if (fDebug) AliInfo("In Forwards terminate");
  TStopwatch swt;
  swt.Start();

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

  // Get our histograms from the container 
  TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
  TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
  TH1I* hEventsAcc   = 0;
  TH1I* hTriggers    = 0;
  if (!fEventInspector.FetchHistograms(list, 
				       hEventsTr, 
				       hEventsTrVtx, 
				       hEventsAcc,
				       hTriggers)) { 
    AliError(Form("Didn't get histograms from event selector "
		  "(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p)", 
		  hEventsTr, hEventsTrVtx,hEventsAcc));
    return;
  }

  TStopwatch swf;
  swf.Start();
  fEnergyFitter.Fit(list2);
  swf.Stop();
  AliInfoF("Fitting took %d real-time seconds, and %f CPU seconds", 
	   Int_t(swf.RealTime()), swf.CpuTime());

  fSharingFilter.Terminate(list,list2,Int_t(hEventsTr->Integral()));
  fDensityCalculator.Terminate(list,list2,Int_t(hEventsTrVtx->Integral()));

  if (fDebug) AliInfoF("Posting post processing results to %s", 
		       list2->GetName());
  PostData(2, list2);

  swt.Stop();
  AliInfoF("Terminate took %d real-time seconds, and %f CPU seconds", 
	   Int_t(swt.RealTime()), swt.CpuTime());

  
}

//____________________________________________________________________
void
AliForwardQATask::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);
  GetEnergyFitter()     .Print(option);
  GetSharingFilter()    .Print(option);
  GetDensityCalculator().Print(option);
  gROOT->DecreaseDirLevel();
}

//
// EOF
//
Commit	Line	Data
96624385	1	//
	2	// Calculate the multiplicity in the forward regions event-by-event
	3	//
	4	// Inputs:
	5	// - AliESDEvent
	6	//
	7	// Outputs:
	8	// - AliAODForwardMult
	9	//
	10	// Histograms
	11	//
	12	// Corrections used
	13	//
	14	#include "AliForwardQATask.h"
	15	#include "AliForwardUtil.h"
	16	#include "AliTriggerAnalysis.h"
	17	#include "AliPhysicsSelection.h"
	18	#include "AliLog.h"
	19	#include "AliESDEvent.h"
	20	#include "AliAODHandler.h"
	21	#include "AliMultiplicity.h"
	22	#include "AliInputEventHandler.h"
	23	#include "AliForwardCorrectionManager.h"
	24	#include "AliAnalysisManager.h"
	25	#include "AliAODForwardMult.h"
	26	#include <TH1.h>
	27	#include <TDirectory.h>
	28	#include <TTree.h>
	29	#include <TROOT.h>
e2213ed5	30	#include <TStopwatch.h>
96624385	31
	32	//====================================================================
	33	AliForwardQATask::AliForwardQATask()
	34	: AliAnalysisTaskSE(),
	35	fEnableLowFlux(false),
	36	fFirstEvent(true),
	37	fCorrManager(0),
	38	fESDFMD(),
	39	fHistos(),
	40	fEventInspector(),
	41	fEnergyFitter(),
	42	fSharingFilter(),
	43	fDensityCalculator(),
9b2f2e39	44	fList(0),
9b2f2e39	45	fDebug(0)
96624385	46	{
	47	//
	48	// Constructor
	49	//
	50	}
	51
	52	//____________________________________________________________________
	53	AliForwardQATask::AliForwardQATask(const char* name)
	54	: AliAnalysisTaskSE(name),
	55	fEnableLowFlux(false),
	56	fFirstEvent(true),
	57	fCorrManager(0),
	58	fESDFMD(),
	59	fHistos(),
	60	fEventInspector("event"),
	61	fEnergyFitter("energy"),
	62	fSharingFilter("sharing"),
	63	fDensityCalculator("density"),
9b2f2e39	64	fList(0),
9b2f2e39	65	fDebug(0)
96624385	66	{
	67	//
	68	// Constructor
	69	//
	70	// Parameters:
	71	// name Name of task
	72	//
	73	DefineOutput(1, TList::Class());
	74	DefineOutput(2, TList::Class());
	75	fCorrManager = &AliForwardCorrectionManager::Instance();
	76	fEnergyFitter.SetNParticles(1); // Just find the 1st peak
	77	fEnergyFitter.SetDoMakeObject(false);
	78	fEnergyFitter.SetUseIncreasingBins(true);
	79	fEnergyFitter.SetDoFits(kTRUE);
	80	fEnergyFitter.SetLowCut(0.4);
	81	fEnergyFitter.SetFitRangeBinWidth(4);
	82	fEnergyFitter.SetMinEntries(1000);
	83	}
	84
	85	//____________________________________________________________________
	86	AliForwardQATask::AliForwardQATask(const AliForwardQATask& o)
	87	: AliAnalysisTaskSE(o),
	88	fEnableLowFlux(o.fEnableLowFlux),
	89	fFirstEvent(o.fFirstEvent),
	90	fCorrManager(o.fCorrManager),
	91	fESDFMD(o.fESDFMD),
	92	fHistos(o.fHistos),
	93	fEventInspector(o.fEventInspector),
	94	fEnergyFitter(o.fEnergyFitter),
	95	fSharingFilter(o.fSharingFilter),
	96	fDensityCalculator(o.fDensityCalculator),
9b2f2e39	97	fList(o.fList),
9b2f2e39	98	fDebug(o.fDebug)
96624385	99	{
	100	//
	101	// Copy constructor
	102	//
	103	// Parameters:
	104	// o Object to copy from
	105	//
	106	DefineOutput(1, TList::Class());
	107	DefineOutput(2, TList::Class());
	108	}
	109
	110	//____________________________________________________________________
	111	AliForwardQATask&
	112	AliForwardQATask::operator=(const AliForwardQATask& o)
	113	{
	114	//
	115	// Assignment operator
	116	//
	117	// Parameters:
	118	// o Object to assign from
	119	//
	120	// Return:
	121	// Reference to this object
	122	//
d015ecfe	123	if (&o == this) return *this;
96624385	124	AliAnalysisTaskSE::operator=(o);
	125
	126	fEnableLowFlux = o.fEnableLowFlux;
	127	fFirstEvent = o.fFirstEvent;
	128	fCorrManager = o.fCorrManager;
	129	fEventInspector = o.fEventInspector;
	130	fEnergyFitter = o.fEnergyFitter;
	131	fSharingFilter = o.fSharingFilter;
	132	fDensityCalculator = o.fDensityCalculator;
	133	fHistos = o.fHistos;
	134	fList = o.fList;
9b2f2e39	135	fDebug = o.fDebug;
96624385	136
	137	return *this;
	138	}
	139
	140	//____________________________________________________________________
	141	void
	142	AliForwardQATask::SetDebug(Int_t dbg)
	143	{
	144	//
	145	// Set debug level
	146	//
	147	// Parameters:
	148	// dbg Debug level
	149	//
9b2f2e39	150	fDebug = dbg;
96624385	151	fEventInspector.SetDebug(dbg);
	152	fEnergyFitter.SetDebug(dbg);
	153	fSharingFilter.SetDebug(dbg);
	154	fDensityCalculator.SetDebug(dbg);
	155	}
	156
	157	//____________________________________________________________________
	158	Bool_t
57522224	159	AliForwardQATask::CheckCorrections(UInt_t what)
96624385	160	{
	161	//
	162	// Check if all needed corrections are there and accounted for. If not,
	163	// do a Fatal exit
	164	//
	165	// Parameters:
	166	// what Which corrections is needed
	167	//
	168	// Return:
	169	// true if all present, false otherwise
	170	//
	171
	172	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	173	// Check that we have the energy loss fits, needed by
	174	// AliFMDSharingFilter
	175	// AliFMDDensityCalculator
8449e3e0	176	if (what & AliForwardCorrectionManager::kELossFits) {
	177	if (!fcm.GetELossFit()) {
	178	AliWarning("No energy loss fits");
	179
	180	// Fall-back values if we do not have the energy loss fits
	181	AliFMDMultCuts& sfLCuts = GetSharingFilter().GetLCuts();
	182	if (sfLCuts.GetMethod() != AliFMDMultCuts::kFixed) {
	183	Double_t cut = 0.3;
	184	AliWarningF("Using fixed cut @ %f for the lower bound "
	185	"of the sharing filter", cut);
	186	sfLCuts.SetMultCuts(cut);
	187	}
	188	AliFMDMultCuts& sfHCuts = GetSharingFilter().GetHCuts();
	189	if (sfHCuts.GetMethod() != AliFMDMultCuts::kFixed) {
	190	Double_t cut = 100;
	191	AliWarningF("Using fixed cut @ %f for the upper bound "
	192	"of the sharing filter", cut);
	193	sfHCuts.SetMultCuts(cut);
	194	}
	195	AliFMDMultCuts& dcCuts = GetDensityCalculator().GetCuts();
	196	if (dcCuts.GetMethod() != AliFMDMultCuts::kFixed) {
	197	Double_t cut = 0.3;
	198	AliWarningF("Using fixed cut @ %f for the lower bound "
	199	"of the density calculator", cut);
	200	dcCuts.SetMultCuts(cut);
	201	}
57522224	202	}
8449e3e0	203	else
8449e3e0	204	fcm.GetELossFit()->CacheBins(GetDensityCalculator().GetMinQuality());
96624385	205	}
8449e3e0	206
96624385	207	return true;
	208	}
	209
	210	//____________________________________________________________________
	211	Bool_t
	212	AliForwardQATask::ReadCorrections(const TAxis*& pe,
	213	const TAxis*& pv,
8449e3e0	214	Bool_t mc,
8449e3e0	215	Bool_t sat)
96624385	216	{
	217	//
	218	// Read corrections
	219	//
	220	//
	221	UInt_t what = AliForwardCorrectionManager::kAll;
	222	what ^= AliForwardCorrectionManager::kDoubleHit;
	223	what ^= AliForwardCorrectionManager::kVertexBias;
	224	what ^= AliForwardCorrectionManager::kAcceptance;
	225	what ^= AliForwardCorrectionManager::kMergingEfficiency;
	226
	227	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
2a50d35b	228	// We allow fall-back queries so that we may proceed in case we have no
	229	// valid corrections
	230	fcm.SetEnableFallBack(true);
8449e3e0	231	if (!fcm.Init(GetEventInspector().GetRunNumber(),
8449e3e0	232	GetEventInspector().GetCollisionSystem(),
96624385	233	GetEventInspector().GetEnergy(),
	234	GetEventInspector().GetField(),
	235	mc,
8449e3e0	236	sat,
	237	what,
	238	false)) return false;
460a5c02	239	if (!CheckCorrections(what)) {
	240	return false;
	241	}
96624385	242
	243	// Sett our persistency pointer
	244	// fCorrManager = &fcm;
	245
	246	// Get the eta axis from the secondary maps - if read in
	247	if (!pe) {
	248	pe = fcm.GetEtaAxis();
	249	if (!pe) AliFatal("No eta axis defined");
	250	}
	251	// Get the vertex axis from the secondary maps - if read in
	252	if (!pv) {
	253	pv = fcm.GetVertexAxis();
	254	if (!pv) AliFatal("No vertex axis defined");
	255	}
	256
	257	return true;
	258	}
	259
	260	//____________________________________________________________________
	261	AliESDEvent*
	262	AliForwardQATask::GetESDEvent()
	263	{
	264	//
	265	// Get the ESD event. IF this is the first event, initialise
	266	//
241cca4d	267	DGUARD(fDebug,2,"Get the ESD event");
	268	if (IsZombie()) {
	269	DMSG(fDebug,3,"We're a Zombie - bailing out");
	270	return 0;
	271	}
96624385	272	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	273	if (!esd) {
	274	AliWarning("No ESD event found for input event");
	275	return 0;
	276	}
	277
	278	// On the first event, initialize the parameters
	279	if (fFirstEvent && esd->GetESDRun()) {
	280	GetEventInspector().ReadRunDetails(esd);
	281
422a78c8	282	AliInfoF("Initializing with parameters from the ESD:\n"
	283	" AliESDEvent::GetBeamEnergy() ->%f\n"
	284	" AliESDEvent::GetBeamType() ->%s\n"
	285	" AliESDEvent::GetCurrentL3() ->%f\n"
	286	" AliESDEvent::GetMagneticField()->%f\n"
6ff251d8	287	" AliESDEvent::GetRunNumber() ->%d",
422a78c8	288	esd->GetBeamEnergy(),
	289	esd->GetBeamType(),
	290	esd->GetCurrentL3(),
	291	esd->GetMagneticField(),
	292	esd->GetRunNumber());
96624385	293
96624385	294
5934a3e3	295	if (!SetupForData()) {
422a78c8	296	AliWarning("Initialisation of sub algorithms failed!");
6ff251d8	297	SetZombie(true);
6ff251d8	298	esd = 0;
422a78c8	299	return 0;
422a78c8	300	}
460a5c02	301	AliInfoF("Clearing first event flag from %s to false",
	302	fFirstEvent ? "true" : "false");
	303	fFirstEvent = false;
96624385	304	}
	305	return esd;
	306	}
	307	//____________________________________________________________________
422a78c8	308	Bool_t
5934a3e3	309	AliForwardQATask::SetupForData()
96624385	310	{
	311	//
	312	// Initialise the sub objects and stuff. Called on first event
	313	//
	314	//
	315	const TAxis* pe = 0;
	316	const TAxis* pv = 0;
	317
d70d2a44	318	if (!ReadCorrections(pe,pv)) {
	319	AliWarning("Using default binning");
	320	pv = new TAxis(10,-10, 10);
	321	pe = new TAxis(240,-6,6);
460a5c02	322	}
96624385	323
	324	fHistos.Init(*pe);
	325
5934a3e3	326	fEventInspector.SetupForData(*pv);
	327	fEnergyFitter.SetupForData(*pe);
	328	fSharingFilter.SetupForData(*pe);
	329	fDensityCalculator.SetupForData(*pe);
96624385	330
96624385	331	this->Print();
422a78c8	332
422a78c8	333	return true;
96624385	334	}
	335
	336	//____________________________________________________________________
	337	void
	338	AliForwardQATask::UserCreateOutputObjects()
	339	{
	340	//
	341	// Create output objects
	342	//
	343	//
	344	fList = new TList;
	345	fList->SetOwner();
	346
5934a3e3	347	fEventInspector.CreateOutputObjects(fList);
	348	fEnergyFitter.CreateOutputObjects(fList);
	349	fSharingFilter.CreateOutputObjects(fList);
	350	fDensityCalculator.CreateOutputObjects(fList);
96624385	351
	352	PostData(1, fList);
	353	}
	354	//____________________________________________________________________
	355	void
	356	AliForwardQATask::UserExec(Option_t*)
	357	{
	358	//
	359	// Process each event
	360	//
	361	// Parameters:
	362	// option Not used
	363	//
241cca4d	364	DGUARD(fDebug,1,"Process the input event");
96624385	365
	366	// static Int_t cnt = 0;
	367	// cnt++;
	368	// Get the input data
	369	AliESDEvent* esd = GetESDEvent();
422a78c8	370	if (!esd) {
	371	AliWarning("Got no ESD event");
	372	return;
	373	}
460a5c02	374	if (fFirstEvent) {
	375	// If the first event flag wasn't cleared in the above call to
	376	// GetESDEvent, we should not do anything, since nothing has been
	377	// initialised yet, so we opt out here (with a warning)
	378	AliWarning("Nothing has been initialized yet, opt'ing out");
	379	return;
	380	}
96624385	381
	382	// Clear stuff
	383	fHistos.Clear();
	384	fESDFMD.Clear();
	385
	386	Bool_t lowFlux = kFALSE;
	387	UInt_t triggers = 0;
	388	UShort_t ivz = 0;
5ca83fee	389	TVector3 ip;
96624385	390	Double_t cent = -1;
	391	UShort_t nClusters = 0;
	392	UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
5ca83fee	393	ivz, ip, cent, nClusters);
96624385	394
241cca4d	395	Bool_t ok = true;
	396	if (found & AliFMDEventInspector::kNoEvent) ok = false;
	397	if (found & AliFMDEventInspector::kNoTriggers) ok = false;
	398	if (found & AliFMDEventInspector::kNoSPD) ok = false;
	399	if (found & AliFMDEventInspector::kNoFMD) ok = false;
	400	if (found & AliFMDEventInspector::kNoVertex) ok = false;
	401	if (triggers & AliAODForwardMult::kPileUp) ok = false;
81775aba	402	if (triggers & AliAODForwardMult::kA) ok = false;
	403	if (triggers & AliAODForwardMult::kC) ok = false;
	404	if (triggers & AliAODForwardMult::kE) ok = false;
	405	if (!(triggers & AliAODForwardMult::kOffline)) ok = false;
241cca4d	406	if (found & AliFMDEventInspector::kBadVertex) ok = false;
	407	if (!ok) {
	408	DMSG(fDebug,2,"Event failed selection: %s",
	409	AliFMDEventInspector::CodeString(found));
	410	return;
	411	}
5ca83fee	412	DMSG(fDebug,2,"Event triggers: %s",
5ca83fee	413	AliAODForwardMult::GetTriggerString(triggers));
96624385	414
	415	// We we do not want to use low flux specific code, we disable it here.
	416	if (!fEnableLowFlux) lowFlux = false;
	417
	418	// Get FMD data
	419	AliESDFMD* esdFMD = esd->GetFMDData();
	420
	421	// Run the energy loss fitter
	422	if (!fEnergyFitter.Accumulate(*esdFMD, cent,
	423	triggers & AliAODForwardMult::kEmpty)) {
	424	AliWarning("Energy fitter failed");
	425	return;
	426	}
4077f3e8	427
96624385	428	// // Apply the sharing filter (or hit merging or clustering if you like)
5ca83fee	429	if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, ip.Z())) {
96624385	430	AliWarning("Sharing filter failed!");
	431	return;
	432	}
4077f3e8	433
96624385	434	// Calculate the inclusive charged particle density
57522224	435	if (!fDensityCalculator.Calculate(fESDFMD, fHistos, lowFlux, cent, ip)) {
96624385	436	// if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) {
	437	AliWarning("Density calculator failed!");
	438	return;
	439	}
4077f3e8	440
96624385	441	PostData(1, fList);
	442	}
	443
	444	//____________________________________________________________________
	445	void
	446	AliForwardQATask::Terminate(Option_t*)
	447	{
	448	//
	449	// End of job
	450	//
	451	// Parameters:
	452	// option Not used
	453	//
9b2f2e39	454	if (fDebug) AliInfo("In Forwards terminate");
e2213ed5	455	TStopwatch swt;
e2213ed5	456	swt.Start();
96624385	457
	458	TList* list = dynamic_cast<TList*>(GetOutputData(1));
	459	if (!list) {
	460	AliError(Form("No output list defined (%p)", GetOutputData(1)));
	461	if (GetOutputData(1)) GetOutputData(1)->Print();
	462	return;
	463	}
	464
5934a3e3	465	// Make a deep copy and post that as output 2
	466	TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults",
	467	list->GetName())));
	468	list2->SetOwner();
	469
96624385	470	// Get our histograms from the container
	471	TH1I* hEventsTr = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
	472	TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
5ca83fee	473	TH1I* hEventsAcc = 0;
96624385	474	TH1I* hTriggers = 0;
5ca83fee	475	if (!fEventInspector.FetchHistograms(list,
	476	hEventsTr,
	477	hEventsTrVtx,
	478	hEventsAcc,
	479	hTriggers)) {
96624385	480	AliError(Form("Didn't get histograms from event selector "
5ca83fee	481	"(hEventsTr=%p,hEventsTrVtx=%p,hEventsAcc=%p)",
5ca83fee	482	hEventsTr, hEventsTrVtx,hEventsAcc));
96624385	483	return;
	484	}
	485
e2213ed5	486	TStopwatch swf;
e2213ed5	487	swf.Start();
57522224	488	fEnergyFitter.Fit(list2);
e2213ed5	489	swf.Stop();
	490	AliInfoF("Fitting took %d real-time seconds, and %f CPU seconds",
	491	Int_t(swf.RealTime()), swf.CpuTime());
96624385	492
5934a3e3	493	fSharingFilter.Terminate(list,list2,Int_t(hEventsTr->Integral()));
5934a3e3	494	fDensityCalculator.Terminate(list,list2,Int_t(hEventsTrVtx->Integral()));
96624385	495
9b2f2e39	496	if (fDebug) AliInfoF("Posting post processing results to %s",
9b2f2e39	497	list2->GetName());
96624385	498	PostData(2, list2);
96624385	499
e2213ed5	500	swt.Stop();
	501	AliInfoF("Terminate took %d real-time seconds, and %f CPU seconds",
	502	Int_t(swt.RealTime()), swt.CpuTime());
	503
5934a3e3	504
96624385	505	}
	506
	507	//____________________________________________________________________
	508	void
	509	AliForwardQATask::Print(Option_t* option) const
	510	{
	511	//
	512	// Print information
	513	//
	514	// Parameters:
	515	// option Not used
	516	//
	517
	518	std::cout << ClassName() << ": " << GetName() << "\n"
	519	<< " Enable low flux code: " << (fEnableLowFlux ? "yes" : "no")
	520	<< "\n"
	521	<< " Off-line trigger mask: 0x"
	522	<< std::hex << std::setfill('0')
	523	<< std::setw (8) << fOfflineTriggerMask
	524	<< std::dec << std::setfill (' ') << std::endl;
	525	gROOT->IncreaseDirLevel();
	526	if (fCorrManager) fCorrManager->Print();
	527	else
	528	std::cout << " Correction manager not set yet" << std::endl;
	529	GetEventInspector() .Print(option);
	530	GetEnergyFitter() .Print(option);
	531	GetSharingFilter() .Print(option);
57522224	532	GetDensityCalculator().Print(option);
96624385	533	gROOT->DecreaseDirLevel();
	534	}
	535
	536	//
	537	// EOF
	538	//