[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskPtEMCalTrigger.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
 * Analysis task of the pt analysis on EMCal-triggered events
 *
 *   Author: Markus Fasel
 */

#include <map>
#include <cstring>
#include <iostream>
#include <memory>
#include <vector>
#include <string>
#include <sstream>

#include <TDirectory.h>
#include <TH1.h>
#include <THashList.h>
#include <TList.h>
#include <TKey.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TString.h>

#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"

#include "AliEMCalHistoContainer.h"
#include "AliAnalysisTaskPtEMCalTrigger.h"

ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)

namespace EMCalTriggerPtAnalysis {

	//______________________________________________________________________________
	AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
                		AliAnalysisTaskSE(),
                		fResults(NULL),
                		fHistos(NULL),
                		fListTrackCuts(NULL)
	{
		/*
		 * Dummy constructor, initialising the values with default (NULL) values
		 */
	}

	//______________________________________________________________________________
	AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
                		AliAnalysisTaskSE(name),
                		fResults(NULL),
                		fHistos(NULL),
                		fListTrackCuts(NULL)
	{
		/*
		 * Main constructor, setting default values for eta and zvertex cut
		 */
		DefineOutput(1, TList::Class());

		fListTrackCuts = new TList;
		fListTrackCuts->SetOwner(false);

		// Set default cuts
		fEtaRange.SetLimits(-0.8, 0.8);

	}

	//______________________________________________________________________________
	AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
		/*
		 * Destructor, deleting output
		 */
		//if(fTrackSelection) delete fTrackSelection;
		if(fHistos) delete fHistos;
		if(fListTrackCuts) delete fListTrackCuts;
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
		/*
		 * Create the list of output objects and define the histograms.
		 * Also adding the track cuts to the list of histograms.
		 */
		fResults = new TList;
		fResults->SetOwner();

		fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
		fHistos->ReleaseOwner();

		std::map<std::string, std::string> triggerCombinations;
		const char *triggernames[6] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal"};
		// Define axes for the trigger correlation histogram
		const TAxis *triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis *) * 5);
		const char *binlabels[2] = {"OFF", "ON"};
		TAxis mytrgaxis[5];
		for(int itrg = 0; itrg < 5; ++itrg){
			DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
			triggeraxis[itrg] = mytrgaxis+itrg;
		}
		// Define names and titles for different triggers in the histogram container
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "jet-triggered events (high threshold)"));
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "jet-triggered events (low threshold)"));
		triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events (low threshold)"));
		// Define axes for the pt histogram
		// Dimensions:
		// 1. pt
		// 2. eta
		// 3. phi
		// 4. vertex
		// 5. pileup (0 = all events, 1 = after pileup rejection)
		// 6. track cuts (0 = no cuts; 1 = after std cuts)
		TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
		pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
		CreateDefaultPtBinning(ptbinning);
		CreateDefaultZVertexBinning(zvertexBinning);
		CreateDefaultEtaBinning(etabinning);
		TAxis htrackaxes[6];
		DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
		DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
		DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
		DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
		DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
		DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
		const TAxis *trackaxes[6];
		for(int iaxis = 0; iaxis < 6; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
		try{
			for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
				const std::string name = it->first, &title = it->second;
				// Create event-based histogram
				fHistos->CreateTH2(Form("hEventHist%s", name.c_str()), Form("Event-based data for %s events; pileup rejection; z_{V} (cm)", title.c_str()), pileupaxis, zvertexBinning);
				// Create track-based histogram
				fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 6, trackaxes);
			}
			fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
		} catch (HistoContainerContentException &e){
			std::stringstream errormessage;
			errormessage << "Creation of histogram failed: " << e.what();
			AliError(errormessage.str().c_str());
		}
		fResults->Add(fHistos->GetListOfHistograms());
		if(fListTrackCuts && fListTrackCuts->GetEntries()){
			TIter cutIter(fListTrackCuts);
			AliESDtrackCuts *cutObject(NULL);
			while((cutObject = dynamic_cast<AliESDtrackCuts *>(cutIter()))){
				cutObject->DefineHistograms();
				fResults->Add(cutObject);
			}
		}
		PostData(1, fResults);
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::UserExec(Option_t* /*option*/){
		/*
		 * Runs the event loop
		 *
		 * @param option: Additional options
		 */

		// Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
		AliESDEvent *esd = static_cast<AliESDEvent *>(fInputEvent);
		const AliESDVertex *vtxTracks = esd->GetPrimaryVertex(),
				*vtxSPD = esd->GetPrimaryVertexSPD();
		if(!(vtxTracks && vtxSPD)) return;
		if(vtxTracks->GetNContributors() < 1 || vtxSPD->GetNContributors() < 1) return;

		double triggers[5]; memset(triggers, 0, sizeof(double) *5);
		if(fInputHandler->IsEventSelected() & AliVEvent::kINT7)
			triggers[0] = 1.;

		std::vector<std::string> triggerstrings;
		if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
			// EMCal-triggered event, distinguish types
			TString trgstr(fInputEvent->GetFiredTriggerClasses());
			if(trgstr.Contains("EJ1")){
				triggerstrings.push_back("EMCJHigh");
				triggers[1] = 1;
			}
			if(trgstr.Contains("EJ2")){
				triggerstrings.push_back("EMCJLow");
				triggers[2] = 1;
			}
			if(trgstr.Contains("EG1")){
				triggerstrings.push_back("EMCGHigh");
				triggers[3] = 1;
			}
			if(trgstr.Contains("EG2")){
				triggerstrings.push_back("EMCGLow");
				triggers[4] = 1;
			}
		}
		try{
			fHistos->FillTHnSparse("hEventTriggers", triggers);
		} catch (HistoContainerContentException &e){
			std::stringstream errormessage;
			errormessage << "Filling of histogram failed: " << e.what();
			AliError(errormessage.str().c_str());
		}

		// apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
		bool isPileupEvent = esd->IsPileupFromSPD();
		isPileupEvent = isPileupEvent || (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
		double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
		isPileupEvent = isPileupEvent || (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);

		// Fill event-based histogram
		const double &zv = vtxTracks->GetZ();
		if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
		if(!triggerstrings.size()) // Non-EMCal-triggered
			FillEventHist("NoEMCal", zv, isPileupEvent);
		else{
			// EMCal-triggered events
			for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
				FillEventHist(it->c_str(), zv, isPileupEvent);
		}

		AliESDtrack *track(NULL);
		// Loop over all tracks (No cuts applied)
		for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
			track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrk));
			// first fill without pielup cut
			if(fEtaRange.IsInRange(track->Eta())) continue;
			if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0);
			if(!triggerstrings.size()) // Non-EMCal-triggered
				FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0);
			else {
				// EMCal-triggered events
				for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
					FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0);
			}
		}

		// Now apply track selection cuts
		// allow for several track selections to be tested at the same time
		// each track selection gets a different cut ID starting from 1
		// cut ID 0 is reserved for the case of no cuts
		if(fListTrackCuts && fListTrackCuts->GetEntries()){
			for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
				AliESDtrackCuts *trackSelection = static_cast<AliESDtrackCuts *>(fListTrackCuts->At(icut));
				std::auto_ptr<TObjArray> acceptedTracks(trackSelection->GetAcceptedTracks(esd));
				TIter trackIter(acceptedTracks.get());
				while((track = dynamic_cast<AliESDtrack *>(trackIter()))){
					if(!fEtaRange.IsInRange(track->Eta())) continue;
					if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1);
					if(!triggerstrings.size()) // Non-EMCal-triggered
						FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1);
					else {
						// EMCal-triggered events
						for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
							FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1);
					}
				}
			}
		}
		PostData(1, fResults);
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
		/*
		 * Creating the default pt binning.
		 *
		 * @param binning: Array where to store the results.
		 */
		std::vector<double> mybinning;
		std::map<double,double> definitions;
		definitions.insert(std::pair<double,double>(2.5, 0.1));
		definitions.insert(std::pair<double,double>(7., 0.25));
		definitions.insert(std::pair<double,double>(15., 0.5));
		definitions.insert(std::pair<double,double>(25., 1.));
		definitions.insert(std::pair<double,double>(40., 2.5));
		definitions.insert(std::pair<double,double>(60., 5.));
		definitions.insert(std::pair<double,double>(100., 5.));
		double currentval = 0;
		for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
			double limit = id->first, binwidth = id->second;
			while(currentval <= limit){
				currentval += binwidth;
				mybinning.push_back(currentval);
			}
		}
		binning.Set(mybinning.size());
		int ib = 0;
		for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
			binning[ib++] = *it;
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
		/*
		 * Creating default z-Vertex binning.
		 *
		 * @param binning: Array where to store the results.
		 */
		std::vector<double> mybinning;
		double currentval = -40;
		mybinning.push_back(currentval);
		while(currentval <= 40.){
			currentval += 1.;
			mybinning.push_back(currentval);
		}
		binning.Set(mybinning.size());
		int ib = 0;
		for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
			binning[ib++] = *it;
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
		/*
		 * Creating default z-Vertex binning.
		 *
		 * @param binning: Array where to store the results.
		 */
		std::vector<double> mybinning;
		double currentval = -0.8;
		mybinning.push_back(currentval);
		while(currentval <= 0.8){
			currentval += 0.1;
			mybinning.push_back(currentval);
		}
		binning.Set(mybinning.size());
		int ib = 0;
		for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
			binning[ib++] = *it;
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
			const char* title, const TArrayD& binning, const char** labels) {
		/*
		 * Define an axis with a given binning
		 *
		 * @param axis: Axis to be defined
		 * @param name: Name of the axis
		 * @param title: Title of the axis
		 * @param binning: axis binning
		 * @param labels (@optional): array of bin labels
		 */
		axis.Set(binning.GetSize()-1, binning.GetArray());
		axis.SetName(name);
		axis.SetTitle(title);
		if(labels){
			for(int ib = 1; ib <= axis.GetNbins(); ++ib)
				axis.SetBinLabel(ib, labels[ib]);
		}
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
			const char* title, int nbins, double min, double max,
			const char** labels) {
		/*
		 * Define an axis with number of bins from min to max
		 *
		 * @param axis: Axis to be defined
		 * @param name: Name of the axis
		 * @param title: Title of the axis
		 * @param nbins: Number of bins
		 * @param min: lower limit of the axis
		 * @param max: upper limit of the axis
		 * @param labels (@optional): array of bin labels
		 */
		axis.Set(nbins, min, max);
		axis.SetName(name);
		axis.SetTitle(title);
		if(labels){
			for(int ib = 1; ib <= axis.GetNbins(); ++ib)
				axis.SetBinLabel(ib, labels[ib]);
		}
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
			double vz, bool isPileup) {
		/*
		 * Fill event-based histogram
		 *
		 * @param trigger: name of the trigger configuration to be processed
		 * @param vz: z-position of the vertex
		 * @param isPileup: signalises if the event is flagged as pileup event
		 */
		char histname[1024];
		sprintf(histname, "hEventHist%s", trigger);
		try{
			fHistos->FillTH2(histname, vz, 0);
		} catch (HistoContainerContentException &e){
			std::stringstream errormessage;
			errormessage << "Filling of histogram failed: " << e.what();
			AliError(errormessage.str().c_str());
		}
		if(!isPileup){
			try{
				fHistos->FillTH2(histname, vz, 1);
			} catch(HistoContainerContentException &e){
				std::stringstream errormessage;
				errormessage << "Filling of histogram failed: " << e.what();
				AliError(errormessage.str().c_str());
			}
		}
	}

	//______________________________________________________________________________
	void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
			const AliESDtrack* track, double vz, bool isPileup, int cut) {
		/*
		 * Fill track-based histogram with corresponding information
		 *
		 * @param trigger: name of the trigger
		 * @param track: ESD track selected
		 * @param vz: z-position of the vertex
		 * @param isPileup: flag event as pileup event
		 * @param cut: id of the cut (0 = no cut)
		 */
		double data[6] = {track->Pt(), track->Eta(), track->Phi(), vz, 0, cut};
		char histname[1024];
		sprintf(histname, "hTrackHist%s", trigger);
		try{
			fHistos->FillTHnSparse(histname, data);
		} catch (HistoContainerContentException &e){
			std::stringstream errormessage;
			errormessage << "Filling of histogram failed: " << e.what();
			AliError(errormessage.str().c_str());
		}
		if(!isPileup){
			data[4] = 1;
			try{
				fHistos->FillTHnSparse(histname, data);
			} catch (HistoContainerContentException &e){
				std::stringstream errormessage;
				errormessage << "Filling of histogram failed: " << e.what();
				AliError(errormessage.str().c_str());
			}
		}
	}

}
Commit	Line	Data
46f589c2	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	* Analysis task of the pt analysis on EMCal-triggered events
	18	*
	19	* Author: Markus Fasel
	20	*/
	21
	22	#include <map>
	23	#include <cstring>
	24	#include <iostream>
	25	#include <memory>
	26	#include <vector>
	27	#include <string>
	28	#include <sstream>
	29
	30	#include <TDirectory.h>
	31	#include <TH1.h>
	32	#include <THashList.h>
eed93770	33	#include <TList.h>
46f589c2	34	#include <TKey.h>
	35	#include <TMath.h>
	36	#include <TObjArray.h>
	37	#include <TString.h>
	38
	39	#include "AliESDEvent.h"
	40	#include "AliESDInputHandler.h"
	41	#include "AliESDtrack.h"
	42	#include "AliESDVertex.h"
	43
	44	#include "AliEMCalHistoContainer.h"
	45	#include "AliAnalysisTaskPtEMCalTrigger.h"
	46
	47	ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskPtEMCalTrigger)
	48
	49	namespace EMCalTriggerPtAnalysis {
	50
	51	//______________________________________________________________________________
	52	AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger():
	53	AliAnalysisTaskSE(),
	54	fResults(NULL),
	55	fHistos(NULL),
	56	fListTrackCuts(NULL)
	57	{
	58	/*
	59	* Dummy constructor, initialising the values with default (NULL) values
	60	*/
	61	}
	62
	63	//______________________________________________________________________________
	64	AliAnalysisTaskPtEMCalTrigger::AliAnalysisTaskPtEMCalTrigger(const char *name):
	65	AliAnalysisTaskSE(name),
	66	fResults(NULL),
	67	fHistos(NULL),
	68	fListTrackCuts(NULL)
	69	{
	70	/*
	71	* Main constructor, setting default values for eta and zvertex cut
	72	*/
	73	DefineOutput(1, TList::Class());
	74
	75	fListTrackCuts = new TList;
	76	fListTrackCuts->SetOwner(false);
	77
	78	// Set default cuts
	79	fEtaRange.SetLimits(-0.8, 0.8);
	80
	81	}
	82
	83	//______________________________________________________________________________
	84	AliAnalysisTaskPtEMCalTrigger::~AliAnalysisTaskPtEMCalTrigger(){
	85	/*
	86	* Destructor, deleting output
	87	*/
	88	//if(fTrackSelection) delete fTrackSelection;
	89	if(fHistos) delete fHistos;
	90	if(fListTrackCuts) delete fListTrackCuts;
	91	}
	92
	93	//______________________________________________________________________________
	94	void AliAnalysisTaskPtEMCalTrigger::UserCreateOutputObjects(){
	95	/*
	96	* Create the list of output objects and define the histograms.
	97	* Also adding the track cuts to the list of histograms.
98	*/
99	fResults = new TList;
100	fResults->SetOwner();
101
102	fHistos = new AliEMCalHistoContainer("PtEMCalTriggerHistograms");
103	fHistos->ReleaseOwner();
104
105	std::map<std::string, std::string> triggerCombinations;
106	const char *triggernames[6] = {"MinBias", "EMCJHigh", "EMCJLow", "EMCGHigh", "EMCGLow", "NoEMCal"};
107	// Define axes for the trigger correlation histogram
108	const TAxis triggeraxis[5]; memset(triggeraxis, 0, sizeof(const TAxis ) * 5);
109	const char *binlabels[2] = {"OFF", "ON"};
110	TAxis mytrgaxis[5];
111	for(int itrg = 0; itrg < 5; ++itrg){
112	DefineAxis(mytrgaxis[itrg], triggernames[itrg], triggernames[itrg], 2, -0.5, 1.5, binlabels);
113	triggeraxis[itrg] = mytrgaxis+itrg;
114	}
115	// Define names and titles for different triggers in the histogram container
116	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[0], "min. bias events"));
117	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[1], "jet-triggered events (high threshold)"));
118	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[2], "jet-triggered events (low threshold)"));
119	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[3], "jet-triggered events (high threshold)"));
120	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[4], "jet-triggered events (low threshold)"));
121	triggerCombinations.insert(std::pair<std::string,std::string>(triggernames[5], "non-EMCal-triggered events (low threshold)"));
122	// Define axes for the pt histogram
123	// Dimensions:
124	// 1. pt
125	// 2. eta
126	// 3. phi
127	// 4. vertex
128	// 5. pileup (0 = all events, 1 = after pileup rejection)
129	// 6. track cuts (0 = no cuts; 1 = after std cuts)
130	TArrayD ptbinning, zvertexBinning, etabinning, pileupaxis(3);
131	pileupaxis[0] = -0.5; pileupaxis[1] = 0.5; pileupaxis[2] = 1.5;
132	CreateDefaultPtBinning(ptbinning);
133	CreateDefaultZVertexBinning(zvertexBinning);
134	CreateDefaultEtaBinning(etabinning);
135	TAxis htrackaxes[6];
136	DefineAxis(htrackaxes[0], "pt", "p_{t} (GeV/c)", ptbinning);
137	DefineAxis(htrackaxes[1], "eta", "#eta", etabinning);
eed93770	138	DefineAxis(htrackaxes[2], "phi", "#phi", 20, 0, 2 * TMath::Pi());
46f589c2	139	DefineAxis(htrackaxes[3], "zvertex", "z_{V} (cm)", zvertexBinning);
	140	DefineAxis(htrackaxes[4], "pileup", "Pileup rejection", 2, -0.5, 1.5);
	141	DefineAxis(htrackaxes[5], "trackcuts", "Track Cuts", (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 1, -0.5, (fListTrackCuts ? fListTrackCuts->GetEntries() : 0) + 0.5);
	142	const TAxis *trackaxes[6];
	143	for(int iaxis = 0; iaxis < 6; ++iaxis) trackaxes[iaxis] = htrackaxes + iaxis;
	144	try{
	145	for(std::map<std::string,std::string>::iterator it = triggerCombinations.begin(); it != triggerCombinations.end(); ++it){
	146	const std::string name = it->first, &title = it->second;
	147	// Create event-based histogram
	148	fHistos->CreateTH2(Form("hEventHist%s", name.c_str()), Form("Event-based data for %s events; pileup rejection; z_{V} (cm)", title.c_str()), pileupaxis, zvertexBinning);
	149	// Create track-based histogram
	150	fHistos->CreateTHnSparse(Form("hTrackHist%s", name.c_str()), Form("Track-based data for %s events", title.c_str()), 6, trackaxes);
	151	}
	152	fHistos->CreateTHnSparse("hEventTriggers", "Trigger type per event", 5, triggeraxis);
	153	} catch (HistoContainerContentException &e){
	154	std::stringstream errormessage;
	155	errormessage << "Creation of histogram failed: " << e.what();
	156	AliError(errormessage.str().c_str());
	157	}
	158	fResults->Add(fHistos->GetListOfHistograms());
	159	if(fListTrackCuts && fListTrackCuts->GetEntries()){
	160	TIter cutIter(fListTrackCuts);
	161	AliESDtrackCuts *cutObject(NULL);
	162	while((cutObject = dynamic_cast<AliESDtrackCuts *>(cutIter()))){
	163	cutObject->DefineHistograms();
	164	fResults->Add(cutObject);
	165	}
	166	}
	167	PostData(1, fResults);
	168	}
	169
	170	//______________________________________________________________________________
	171	void AliAnalysisTaskPtEMCalTrigger::UserExec(Option_t* /option/){
	172	/*
	173	* Runs the event loop
	174	*
	175	* @param option: Additional options
	176	*/
	177
	178	// Common checks: Have SPD vertex and primary vertex from tracks, and both need to have at least one contributor
	179	AliESDEvent esd = static_cast<AliESDEvent >(fInputEvent);
	180	const AliESDVertex *vtxTracks = esd->GetPrimaryVertex(),
	181	*vtxSPD = esd->GetPrimaryVertexSPD();
	182	if(!(vtxTracks && vtxSPD)) return;
	183	if(vtxTracks->GetNContributors() < 1 \|\| vtxSPD->GetNContributors() < 1) return;
	184
	185	double triggers[5]; memset(triggers, 0, sizeof(double) *5);
	186	if(fInputHandler->IsEventSelected() & AliVEvent::kINT7)
	187	triggers[0] = 1.;
	188
	189	std::vector<std::string> triggerstrings;
	190	if(fInputHandler->IsEventSelected() & AliVEvent::kEMC7){
	191	// EMCal-triggered event, distinguish types
	192	TString trgstr(fInputEvent->GetFiredTriggerClasses());
	193	if(trgstr.Contains("EJ1")){
	194	triggerstrings.push_back("EMCJHigh");
	195	triggers[1] = 1;
	196	}
	197	if(trgstr.Contains("EJ2")){
	198	triggerstrings.push_back("EMCJLow");
	199	triggers[2] = 1;
	200	}
	201	if(trgstr.Contains("EG1")){
	202	triggerstrings.push_back("EMCGHigh");
203	triggers[3] = 1;
204	}
205	if(trgstr.Contains("EG2")){
206	triggerstrings.push_back("EMCGLow");
207	triggers[4] = 1;
208	}
209	}
210	try{
211	fHistos->FillTHnSparse("hEventTriggers", triggers);
212	} catch (HistoContainerContentException &e){
213	std::stringstream errormessage;
214	errormessage << "Filling of histogram failed: " << e.what();
215	AliError(errormessage.str().c_str());
216	}
217
218	// apply event selection: Combine the Pileup cut from SPD with the other pA Vertex selection cuts.
219	bool isPileupEvent = esd->IsPileupFromSPD();
220	isPileupEvent = isPileupEvent \|\| (TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) > 0.5);
221	double covSPD[6]; vtxSPD->GetCovarianceMatrix(covSPD);
222	isPileupEvent = isPileupEvent \|\| (TString(vtxSPD->GetTitle()).Contains("vertexer:Z") && TMath::Sqrt(covSPD[5]) > 0.25);
223
224	// Fill event-based histogram
225	const double &zv = vtxTracks->GetZ();
226	if(triggers[0]) FillEventHist("MinBias", zv, isPileupEvent);
227	if(!triggerstrings.size()) // Non-EMCal-triggered
228	FillEventHist("NoEMCal", zv, isPileupEvent);
229	else{
230	// EMCal-triggered events
231	for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
232	FillEventHist(it->c_str(), zv, isPileupEvent);
233	}
234
235	AliESDtrack *track(NULL);
236	// Loop over all tracks (No cuts applied)
237	for(int itrk = 0; itrk < fInputEvent->GetNumberOfTracks(); ++itrk){
238	track = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(itrk));
239	// first fill without pielup cut
240	if(fEtaRange.IsInRange(track->Eta())) continue;
241	if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, 0);
242	if(!triggerstrings.size()) // Non-EMCal-triggered
243	FillTrackHist("NoEMCal", track, zv, isPileupEvent, 0);
244	else {
245	// EMCal-triggered events
246	for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
247	FillTrackHist(it->c_str(), track, zv, isPileupEvent, 0);
248	}
249	}
250
251	// Now apply track selection cuts
252	// allow for several track selections to be tested at the same time
253	// each track selection gets a different cut ID starting from 1
254	// cut ID 0 is reserved for the case of no cuts
255	if(fListTrackCuts && fListTrackCuts->GetEntries()){
256	for(int icut = 0; icut < fListTrackCuts->GetEntries(); icut++){
257	AliESDtrackCuts trackSelection = static_cast<AliESDtrackCuts >(fListTrackCuts->At(icut));
258	std::auto_ptr<TObjArray> acceptedTracks(trackSelection->GetAcceptedTracks(esd));
259	TIter trackIter(acceptedTracks.get());
260	while((track = dynamic_cast<AliESDtrack *>(trackIter()))){
261	if(!fEtaRange.IsInRange(track->Eta())) continue;
262	if(triggers[0]) FillTrackHist("MinBias", track, zv, isPileupEvent, icut + 1);
263	if(!triggerstrings.size()) // Non-EMCal-triggered
264	FillTrackHist("NoEMCal", track, zv, isPileupEvent, icut + 1);
265	else {
266	// EMCal-triggered events
267	for(std::vector<std::string>::iterator it = triggerstrings.begin(); it != triggerstrings.end(); ++it)
268	FillTrackHist(it->c_str(), track, zv, isPileupEvent, icut + 1);
269	}
270	}
271	}
272	}
273	PostData(1, fResults);
274	}
275
276	//______________________________________________________________________________
277	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultPtBinning(TArrayD &binning) const{
278	/*
279	* Creating the default pt binning.
280	*
281	* @param binning: Array where to store the results.
282	*/
283	std::vector<double> mybinning;
284	std::map<double,double> definitions;
285	definitions.insert(std::pair<double,double>(2.5, 0.1));
286	definitions.insert(std::pair<double,double>(7., 0.25));
287	definitions.insert(std::pair<double,double>(15., 0.5));
288	definitions.insert(std::pair<double,double>(25., 1.));
289	definitions.insert(std::pair<double,double>(40., 2.5));
290	definitions.insert(std::pair<double,double>(60., 5.));
291	definitions.insert(std::pair<double,double>(100., 5.));
292	double currentval = 0;
293	for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
294	double limit = id->first, binwidth = id->second;
295	while(currentval <= limit){
296	currentval += binwidth;
297	mybinning.push_back(currentval);
298	}
299	}
300	binning.Set(mybinning.size());
301	int ib = 0;
302	for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
303	binning[ib++] = *it;
304	}
305
306	//______________________________________________________________________________
307	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultZVertexBinning(TArrayD &binning) const {
308	/*
309	* Creating default z-Vertex binning.
310	*
311	* @param binning: Array where to store the results.
312	*/
313	std::vector<double> mybinning;
314	double currentval = -40;
315	mybinning.push_back(currentval);
316	while(currentval <= 40.){
eed93770	317	currentval += 1.;
46f589c2	318	mybinning.push_back(currentval);
	319	}
	320	binning.Set(mybinning.size());
	321	int ib = 0;
	322	for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
	323	binning[ib++] = *it;
	324	}
	325
	326	//______________________________________________________________________________
	327	void AliAnalysisTaskPtEMCalTrigger::CreateDefaultEtaBinning(TArrayD& binning) const {
	328	/*
	329	* Creating default z-Vertex binning.
	330	*
	331	* @param binning: Array where to store the results.
	332	*/
	333	std::vector<double> mybinning;
	334	double currentval = -0.8;
	335	mybinning.push_back(currentval);
	336	while(currentval <= 0.8){
eed93770	337	currentval += 0.1;
46f589c2	338	mybinning.push_back(currentval);
	339	}
	340	binning.Set(mybinning.size());
	341	int ib = 0;
	342	for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
	343	binning[ib++] = *it;
	344	}
	345
	346	//______________________________________________________________________________
	347	void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
	348	const char* title, const TArrayD& binning, const char** labels) {
	349	/*
	350	* Define an axis with a given binning
	351	*
	352	* @param axis: Axis to be defined
	353	* @param name: Name of the axis
	354	* @param title: Title of the axis
	355	* @param binning: axis binning
	356	* @param labels (@optional): array of bin labels
	357	*/
	358	axis.Set(binning.GetSize()-1, binning.GetArray());
	359	axis.SetName(name);
	360	axis.SetTitle(title);
	361	if(labels){
	362	for(int ib = 1; ib <= axis.GetNbins(); ++ib)
	363	axis.SetBinLabel(ib, labels[ib]);
	364	}
	365	}
	366
	367	//______________________________________________________________________________
	368	void AliAnalysisTaskPtEMCalTrigger::DefineAxis(TAxis& axis, const char* name,
	369	const char* title, int nbins, double min, double max,
	370	const char** labels) {
	371	/*
	372	* Define an axis with number of bins from min to max
	373	*
	374	* @param axis: Axis to be defined
	375	* @param name: Name of the axis
	376	* @param title: Title of the axis
	377	* @param nbins: Number of bins
	378	* @param min: lower limit of the axis
	379	* @param max: upper limit of the axis
	380	* @param labels (@optional): array of bin labels
	381	*/
	382	axis.Set(nbins, min, max);
	383	axis.SetName(name);
	384	axis.SetTitle(title);
	385	if(labels){
	386	for(int ib = 1; ib <= axis.GetNbins(); ++ib)
	387	axis.SetBinLabel(ib, labels[ib]);
	388	}
	389	}
	390
	391	//______________________________________________________________________________
	392	void AliAnalysisTaskPtEMCalTrigger::FillEventHist(const char* trigger,
	393	double vz, bool isPileup) {
	394	/*
	395	* Fill event-based histogram
	396	*
	397	* @param trigger: name of the trigger configuration to be processed
	398	* @param vz: z-position of the vertex
	399	* @param isPileup: signalises if the event is flagged as pileup event
	400	*/
	401	char histname[1024];
402	sprintf(histname, "hEventHist%s", trigger);
403	try{
404	fHistos->FillTH2(histname, vz, 0);
405	} catch (HistoContainerContentException &e){
406	std::stringstream errormessage;
407	errormessage << "Filling of histogram failed: " << e.what();
408	AliError(errormessage.str().c_str());
409	}
410	if(!isPileup){
411	try{
412	fHistos->FillTH2(histname, vz, 1);
413	} catch(HistoContainerContentException &e){
414	std::stringstream errormessage;
415	errormessage << "Filling of histogram failed: " << e.what();
416	AliError(errormessage.str().c_str());
417	}
418	}
419	}
420
421	//______________________________________________________________________________
422	void AliAnalysisTaskPtEMCalTrigger::FillTrackHist(const char* trigger,
423	const AliESDtrack* track, double vz, bool isPileup, int cut) {
424	/*
425	* Fill track-based histogram with corresponding information
426	*
427	* @param trigger: name of the trigger
428	* @param track: ESD track selected
429	* @param vz: z-position of the vertex
430	* @param isPileup: flag event as pileup event
431	* @param cut: id of the cut (0 = no cut)
432	*/
433	double data[6] = {track->Pt(), track->Eta(), track->Phi(), vz, 0, cut};
434	char histname[1024];
435	sprintf(histname, "hTrackHist%s", trigger);
436	try{
437	fHistos->FillTHnSparse(histname, data);
438	} catch (HistoContainerContentException &e){
439	std::stringstream errormessage;
440	errormessage << "Filling of histogram failed: " << e.what();
441	AliError(errormessage.str().c_str());
442	}
443	if(!isPileup){
444	data[4] = 1;
445	try{
446	fHistos->FillTHnSparse(histname, data);
447	} catch (HistoContainerContentException &e){
448	std::stringstream errormessage;
449	errormessage << "Filling of histogram failed: " << e.what();
450	AliError(errormessage.str().c_str());
451	}
452	}
453	}
454
455	}
456