[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisTaskTotEt.cxx

//_________________________________________________________________________
//  Utility Class for transverse energy studies
//  Task for analysis
//  - reconstruction and MC output
// implementation file
//
//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
//_________________________________________________________________________
//Necessary to read config macros
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>

#include "TChain.h"
#include "TList.h"
#include "TFile.h"
#include "TH2F.h"
#include "THnSparse.h"

#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDtrackCuts.h"

#include "AliAnalysisTaskTotEt.h"
#include "AliAnalysisEtReconstructedPhos.h"
#include "AliAnalysisEtReconstructedEmcal.h"
#include "AliAnalysisEtMonteCarloPhos.h"
#include "AliAnalysisEtMonteCarloEmcal.h"
#include "AliAnalysisEmEtMonteCarlo.h"
#include "AliAnalysisEmEtReconstructed.h"

#include <iostream>
#include <AliCentrality.h>

using namespace std;

ClassImp(AliAnalysisTaskTotEt)

//________________________________________________________________________
AliAnalysisTaskTotEt::AliAnalysisTaskTotEt(const char *name, Bool_t isMc) :
AliAnalysisTaskTransverseEnergy(name, isMc)
,fRecAnalysis(0)
,fMCAnalysis(0)
,fSparseHistRecVsMc(0)
,fSparseRecVsMc(0)
{
    // Constructor
	
    // select if we should use EMCal or PHOS class
    // PHOS by default, EMCal if name string contains EMC
    TString t(name);
    //t.ToUpper();
    if (t.Contains("EMC")) {
		if (t.Contains("Detail")) {
			fMCAnalysis = new AliAnalysisEmEtMonteCarlo();
			fMCAnalysis->SetDataSet(2010);
			fMCAnalysis->Init();
			
			cout << "Instantiating AliAnalysisEmEtMonteCarlo class..."<< endl;
		}
        else if (fMCConfigFile.Length()) {
            cout<<"Rereading AliAnalysisEtMonteCarloEmcal configuration file..."<<endl;
            gROOT->LoadMacro(fMCConfigFile);
            fMCAnalysis = (AliAnalysisEtMonteCarloEmcal *) gInterpreter->ProcessLine("ConfigEtMonteCarlo()");
        }
		
		if (t.Contains("Detail")) {
			fRecAnalysis = new AliAnalysisEmEtReconstructed();
			fRecAnalysis->SetDataSet(2010);
			
			TFile *infile = new TFile("corrections.root");
			fRecAnalysis->SetCorrections((AliAnalysisHadEtCorrections *)infile->Get("hadCorrectionEMCAL"));
			
			fRecAnalysis->Init();
			
			cout << "Instantiating AliAnalysisEmEtReconstructed class..."<< endl;
		}
        else if (fRecoConfigFile.Length()) {
            cout<<"Rereading AliAnalysisEtReconstructedEmcal configuration file..."<<endl;
            gROOT->LoadMacro(fRecoConfigFile);
            fRecAnalysis = (AliAnalysisEtReconstructedEmcal *) gInterpreter->ProcessLine("ConfigEtReconstructed()");
        }
    }
    else {
        if (fMCConfigFile.Length()) {
            cout<<"Rereading AliAnalysisEtMonteCarloPhos configuration file..."<<endl;
            gROOT->LoadMacro(fMCConfigFile);
			
            fMCAnalysis = (AliAnalysisEtMonteCarloPhos *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(false)");
			cout << fMCAnalysis << endl;
        }
		
        if (fRecoConfigFile.Length()) {
            cout<<"Rereading AliAnalysisEtReconstructedPhos configuration file..."<<endl;
            gROOT->LoadMacro(fRecoConfigFile);
            fRecAnalysis = (AliAnalysisEtReconstructedPhos *) gInterpreter->ProcessLine("ConfigEtReconstructed(false)");
        }
    }
    // Define input and output slots here
    // Input slot #0 works with a TChain
    DefineInput(0, TChain::Class());
    // Output slot #1 writes into a TH1 container
	
    DefineOutput(1, TList::Class());
	
}
AliAnalysisTaskTotEt::~AliAnalysisTaskTotEt() {//Destructor
	//    fOutputList->Clear();
    delete fRecAnalysis;
    delete fMCAnalysis;
}

//________________________________________________________________________
void AliAnalysisTaskTotEt::UserCreateOutputObjects()
{
    // Create histograms
    // Called once
	if (fMCAnalysis)
		fMCAnalysis->CreateHistograms();
    fRecAnalysis->CreateHistograms();
    fOutputList = new TList;
    fOutputList->SetOwner();
    fRecAnalysis->FillOutputList(fOutputList);
	if (fMCAnalysis)
		fMCAnalysis->FillOutputList(fOutputList);
    fHistEtRecvsEtMC = new TH2F("fHistEtRecvsEtMC", "Reconstructed E_{T} vs MC E_{T}", 1000, 0.000, 100, 1000, 0.0001, 100);
    fHistEtRecOverEtMC = new TH2F("fHistEtRecOverEtMC", "Reconstructed E_{T} over MC E_{T} vs centrality", 1000, 0.00, 2.0, 11, -0.5, 10.5); 
    fHistDiffEtRecEtMCOverEtMC = new TH2F("fHistDiffEtRecEtMCOverEtMC", "fHistDiffEtRecEtMCOverEtMC", 10000, 0.0, 1000, 1000, -5, 5); 
    fOutputList->Add(fHistEtRecvsEtMC);
    fOutputList->Add(fHistEtRecOverEtMC);
    fOutputList->Add(fHistDiffEtRecEtMCOverEtMC);
	
    Bool_t selectPrimaries=kTRUE;
    if (fRecAnalysis->DataSet()==2009) {
        cout<<"Setting track cuts for the 2009 p+p collisions at 900 GeV"<<endl;
        fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
        fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
        fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
        fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
        //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
        fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
        fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
    }
    if (fRecAnalysis->DataSet()==2010) {
        cout<<"Setting track cuts for the 2010 p+p collisions at 7 GeV"<<endl;
        //cout<<"Warning:  Have not set 2010 track cuts yet!!"<<endl;
        fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
        fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
        fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
        fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
        //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
        fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
        fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
    }
	
    fOutputList->Add(fEsdtrackCutsITSTPC);
    fOutputList->Add(fEsdtrackCutsTPC);
    fOutputList->Add(fEsdtrackCutsITS);
    if (fEsdtrackCutsITSTPC && fEsdtrackCutsTPC) {
        fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
		if (fMCAnalysis)
			fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
        fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
		if (fMCAnalysis)
			fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
        fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
		if (fMCAnalysis)
			fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
        //add ITS stuff!
    }
    else {
        Printf("Error: no track cuts!");
    }
	
}

//________________________________________________________________________
void AliAnalysisTaskTotEt::UserExec(Option_t *)
{ // execute method
	
    fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
    if (!fESDEvent)
    {
        Printf("ERROR: Could not retrieve event");
        return;
    }
	
    //Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber());
	
    AliCentrality *cent = GetCentralityObject();
	
    //if (res == 0 && cent)
    //{
	if (IsPhysicsSelected())
	{
		fRecAnalysis->SetCentralityObject(cent);
		fRecAnalysis->AnalyseEvent(fESDEvent);
		
		AliMCEvent* mcEvent = MCEvent();
		if (mcEvent)
		{
			fMCAnalysis->SetCentralityObject(cent);
			fMCAnalysis->AnalyseEvent(mcEvent, fESDEvent);
			//fMCAnalysis->AnalyseEvent(mcEvent);
		}
		if(fMCAnalysis)
	    {
			fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEtAcc(), fMCAnalysis->GetTotNeutralEtAcc());
			if(fMCAnalysis->GetTotNeutralEtAcc()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEtAcc(), cent->GetCentralityClass10("V0M"));
			if(fMCAnalysis->GetTotNeutralEtAcc()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
	    }
	}
    //}
    // Post output data.
    PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnalysisTaskTotEt::Terminate(Option_t *)
{
    // Draw result to the screen
    // Called once at the end of the query
	
    fOutputList = dynamic_cast<TList*> (GetOutputData(1));
    if (!fOutputList) {
        printf("ERROR: Output list not available\n");
        return;
    }
}
Commit	Line	Data
cf6522d1	1	//_________________________________________________________________________
	2	// Utility Class for transverse energy studies
	3	// Task for analysis
	4	// - reconstruction and MC output
	5	// implementation file
	6	//
	7	//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
	8	//_________________________________________________________________________
964c8159	9	//Necessary to read config macros
	10	#include <TROOT.h>
	11	#include <TSystem.h>
	12	#include <TInterpreter.h>
cf6522d1	13
2fbf38ac	14	#include "TChain.h"
cf6522d1	15	#include "TList.h"
1815c8d7	16	#include "TFile.h"
2fbf38ac	17	#include "TH2F.h"
a403aff5	18	#include "THnSparse.h"
2fbf38ac	19
2fbf38ac	20	#include "AliESDEvent.h"
2fbf38ac	21	#include "AliMCEvent.h"
b5821c13	22	#include "AliESDtrackCuts.h"
cf6522d1	23
2fbf38ac	24	#include "AliAnalysisTaskTotEt.h"
	25	#include "AliAnalysisEtReconstructedPhos.h"
	26	#include "AliAnalysisEtReconstructedEmcal.h"
	27	#include "AliAnalysisEtMonteCarloPhos.h"
	28	#include "AliAnalysisEtMonteCarloEmcal.h"
1815c8d7	29	#include "AliAnalysisEmEtMonteCarlo.h"
1815c8d7	30	#include "AliAnalysisEmEtReconstructed.h"
2fbf38ac	31
2fbf38ac	32	#include <iostream>
a403aff5	33	#include <AliCentrality.h>
2fbf38ac	34
	35	using namespace std;
	36
	37	ClassImp(AliAnalysisTaskTotEt)
	38
	39	//________________________________________________________________________
d3d7bfe9	40	AliAnalysisTaskTotEt::AliAnalysisTaskTotEt(const char *name, Bool_t isMc) :
1815c8d7	41	AliAnalysisTaskTransverseEnergy(name, isMc)
	42	,fRecAnalysis(0)
	43	,fMCAnalysis(0)
	44	,fSparseHistRecVsMc(0)
	45	,fSparseRecVsMc(0)
2fbf38ac	46	{
2fbf38ac	47	// Constructor
1815c8d7	48
2fbf38ac	49	// select if we should use EMCal or PHOS class
	50	// PHOS by default, EMCal if name string contains EMC
	51	TString t(name);
1815c8d7	52	//t.ToUpper();
2fbf38ac	53	if (t.Contains("EMC")) {
1815c8d7	54	if (t.Contains("Detail")) {
	55	fMCAnalysis = new AliAnalysisEmEtMonteCarlo();
	56	fMCAnalysis->SetDataSet(2010);
	57	fMCAnalysis->Init();
	58
	59	cout << "Instantiating AliAnalysisEmEtMonteCarlo class..."<< endl;
	60	}
	61	else if (fMCConfigFile.Length()) {
d3d7bfe9	62	cout<<"Rereading AliAnalysisEtMonteCarloEmcal configuration file..."<<endl;
	63	gROOT->LoadMacro(fMCConfigFile);
	64	fMCAnalysis = (AliAnalysisEtMonteCarloEmcal *) gInterpreter->ProcessLine("ConfigEtMonteCarlo()");
	65	}
1815c8d7	66
	67	if (t.Contains("Detail")) {
	68	fRecAnalysis = new AliAnalysisEmEtReconstructed();
	69	fRecAnalysis->SetDataSet(2010);
	70
	71	TFile *infile = new TFile("corrections.root");
	72	fRecAnalysis->SetCorrections((AliAnalysisHadEtCorrections *)infile->Get("hadCorrectionEMCAL"));
	73
	74	fRecAnalysis->Init();
	75
	76	cout << "Instantiating AliAnalysisEmEtReconstructed class..."<< endl;
	77	}
	78	else if (fRecoConfigFile.Length()) {
d3d7bfe9	79	cout<<"Rereading AliAnalysisEtReconstructedEmcal configuration file..."<<endl;
	80	gROOT->LoadMacro(fRecoConfigFile);
	81	fRecAnalysis = (AliAnalysisEtReconstructedEmcal *) gInterpreter->ProcessLine("ConfigEtReconstructed()");
	82	}
2fbf38ac	83	}
d3d7bfe9	84	else {
	85	if (fMCConfigFile.Length()) {
	86	cout<<"Rereading AliAnalysisEtMonteCarloPhos configuration file..."<<endl;
	87	gROOT->LoadMacro(fMCConfigFile);
1815c8d7	88
a403aff5	89	fMCAnalysis = (AliAnalysisEtMonteCarloPhos *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(false)");
1815c8d7	90	cout << fMCAnalysis << endl;
d3d7bfe9	91	}
1815c8d7	92
d3d7bfe9	93	if (fRecoConfigFile.Length()) {
	94	cout<<"Rereading AliAnalysisEtReconstructedPhos configuration file..."<<endl;
	95	gROOT->LoadMacro(fRecoConfigFile);
a403aff5	96	fRecAnalysis = (AliAnalysisEtReconstructedPhos *) gInterpreter->ProcessLine("ConfigEtReconstructed(false)");
d3d7bfe9	97	}
2fbf38ac	98	}
2fbf38ac	99	// Define input and output slots here
	100	// Input slot #0 works with a TChain
	101	DefineInput(0, TChain::Class());
	102	// Output slot #1 writes into a TH1 container
1815c8d7	103
2fbf38ac	104	DefineOutput(1, TList::Class());
1815c8d7	105
2fbf38ac	106	}
d3d7bfe9	107	AliAnalysisTaskTotEt::~AliAnalysisTaskTotEt() {//Destructor
1815c8d7	108	// fOutputList->Clear();
d3d7bfe9	109	delete fRecAnalysis;
d3d7bfe9	110	delete fMCAnalysis;
951efd81	111	}
2fbf38ac	112
	113	//________________________________________________________________________
	114	void AliAnalysisTaskTotEt::UserCreateOutputObjects()
	115	{
	116	// Create histograms
	117	// Called once
1815c8d7	118	if (fMCAnalysis)
1815c8d7	119	fMCAnalysis->CreateHistograms();
2fbf38ac	120	fRecAnalysis->CreateHistograms();
2fbf38ac	121	fOutputList = new TList;
951efd81	122	fOutputList->SetOwner();
2fbf38ac	123	fRecAnalysis->FillOutputList(fOutputList);
1815c8d7	124	if (fMCAnalysis)
1815c8d7	125	fMCAnalysis->FillOutputList(fOutputList);
a403aff5	126	fHistEtRecvsEtMC = new TH2F("fHistEtRecvsEtMC", "Reconstructed E_{T} vs MC E_{T}", 1000, 0.000, 100, 1000, 0.0001, 100);
	127	fHistEtRecOverEtMC = new TH2F("fHistEtRecOverEtMC", "Reconstructed E_{T} over MC E_{T} vs centrality", 1000, 0.00, 2.0, 11, -0.5, 10.5);
	128	fHistDiffEtRecEtMCOverEtMC = new TH2F("fHistDiffEtRecEtMCOverEtMC", "fHistDiffEtRecEtMCOverEtMC", 10000, 0.0, 1000, 1000, -5, 5);
2fbf38ac	129	fOutputList->Add(fHistEtRecvsEtMC);
a403aff5	130	fOutputList->Add(fHistEtRecOverEtMC);
a403aff5	131	fOutputList->Add(fHistDiffEtRecEtMCOverEtMC);
1815c8d7	132
d3d7bfe9	133	Bool_t selectPrimaries=kTRUE;
	134	if (fRecAnalysis->DataSet()==2009) {
	135	cout<<"Setting track cuts for the 2009 p+p collisions at 900 GeV"<<endl;
	136	fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
	137	fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
	138	fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	139	fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
	140	//ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
	141	fEsdtrackCutsITS = AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
	142	fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
	143	}
	144	if (fRecAnalysis->DataSet()==2010) {
0651f6b4	145	cout<<"Setting track cuts for the 2010 p+p collisions at 7 GeV"<<endl;
0651f6b4	146	//cout<<"Warning: Have not set 2010 track cuts yet!!"<<endl;
d3d7bfe9	147	fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
	148	fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
	149	fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
	150	fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
	151	//ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
0651f6b4	152	fEsdtrackCutsITS = AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
d3d7bfe9	153	fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
d3d7bfe9	154	}
1815c8d7	155
d3d7bfe9	156	fOutputList->Add(fEsdtrackCutsITSTPC);
	157	fOutputList->Add(fEsdtrackCutsTPC);
	158	fOutputList->Add(fEsdtrackCutsITS);
	159	if (fEsdtrackCutsITSTPC && fEsdtrackCutsTPC) {
	160	fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
1815c8d7	161	if (fMCAnalysis)
1815c8d7	162	fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
d3d7bfe9	163	fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
1815c8d7	164	if (fMCAnalysis)
1815c8d7	165	fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
d3d7bfe9	166	fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
1815c8d7	167	if (fMCAnalysis)
1815c8d7	168	fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
d3d7bfe9	169	//add ITS stuff!
	170	}
	171	else {
	172	Printf("Error: no track cuts!");
	173	}
1815c8d7	174
2fbf38ac	175	}
	176
	177	//________________________________________________________________________
	178	void AliAnalysisTaskTotEt::UserExec(Option_t *)
cf6522d1	179	{ // execute method
1815c8d7	180
d3d7bfe9	181	fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
	182	if (!fESDEvent)
	183	{
2fbf38ac	184	Printf("ERROR: Could not retrieve event");
	185	return;
	186	}
1815c8d7	187
	188	//Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber());
	189
b6dd6ad2	190	AliCentrality *cent = GetCentralityObject();
1815c8d7	191
	192	//if (res == 0 && cent)
	193	//{
	194	if (IsPhysicsSelected())
	195	{
	196	fRecAnalysis->SetCentralityObject(cent);
	197	fRecAnalysis->AnalyseEvent(fESDEvent);
	198
	199	AliMCEvent* mcEvent = MCEvent();
	200	if (mcEvent)
	201	{
	202	fMCAnalysis->SetCentralityObject(cent);
	203	fMCAnalysis->AnalyseEvent(mcEvent, fESDEvent);
	204	//fMCAnalysis->AnalyseEvent(mcEvent);
	205	}
	206	if(fMCAnalysis)
a403aff5	207	{
1815c8d7	208	fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEtAcc(), fMCAnalysis->GetTotNeutralEtAcc());
	209	if(fMCAnalysis->GetTotNeutralEtAcc()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEtAcc(), cent->GetCentralityClass10("V0M"));
	210	if(fMCAnalysis->GetTotNeutralEtAcc()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
a403aff5	211	}
1815c8d7	212	}
1815c8d7	213	//}
d3d7bfe9	214	// Post output data.
2fbf38ac	215	PostData(1, fOutputList);
2fbf38ac	216	}
	217
	218	//________________________________________________________________________
	219	void AliAnalysisTaskTotEt::Terminate(Option_t *)
	220	{
	221	// Draw result to the screen
	222	// Called once at the end of the query
1815c8d7	223
2fbf38ac	224	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
	225	if (!fOutputList) {
	226	printf("ERROR: Output list not available\n");
	227	return;
	228	}
	229	}
	230
	231
	232