]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/totEt/AliAnalysisTaskTotEt.cxx
Adding fix for PHOS clusters and adding a couple diagnostic histograms
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisTaskTotEt.cxx
CommitLineData
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)
ee1083d5 8 //_________________________________________________________________________
9 //Necessary to read config macros
964c8159 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
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"
30#include "AliAnalysisEmEtReconstructed.h"
9a365626 31#include "AliAnalysisManager.h"
32#include "AliPIDResponse.h"
33#include "AliTPCPIDResponse.h"
34#include "AliInputEventHandler.h"
2fbf38ac 35
36#include <iostream>
a403aff5 37#include <AliCentrality.h>
2fbf38ac 38
ee1083d5 39 using namespace std;
2fbf38ac 40
41ClassImp(AliAnalysisTaskTotEt)
42
43//________________________________________________________________________
ee1083d5 44 AliAnalysisTaskTotEt::AliAnalysisTaskTotEt(const char *name, Bool_t isMc) :
45 AliAnalysisTaskTransverseEnergy(name, isMc)
9a365626 46 ,fPIDResponse(0)
ee1083d5 47 ,fRecAnalysis(0)
48 ,fMCAnalysis(0)
01b73fb0 49 // ,fSparseHistRecVsMc(0)
50 //,fSparseRecVsMc(0)
2fbf38ac 51{
ee1083d5 52 // Constructor
ee1083d5 53 // select if we should use EMCal or PHOS class
54 // PHOS by default, EMCal if name string contains EMC
55 TString t(name);
c31562f7 56 gROOT->LoadMacro(fMCConfigFile);
57 gROOT->LoadMacro(fRecoConfigFile);
58 //There is a weird problem where the name reverts to a default using the plugin
59 //these lines solve it - there is a function written into ConfigEtMonteCarlo.C which solves this
60 Bool_t isEMCal = t.Contains("EMC");
61 if (!(t.Contains("EMC")) && !(t.Contains("PHOS"))) {//the name does not contain either EMCal or PHOS
62 cout<<"Default arguments called. Reading config file."<<endl;
63 isEMCal = (Bool_t) gInterpreter->ProcessLine("GetIsEMCAL()");
64 isMc = (Bool_t) gInterpreter->ProcessLine("GetIsMC()");
65
66 }
67 //cout<<__FILE__<<" My name is "<<name<<endl;
ee1083d5 68 //t.ToUpper();
c31562f7 69 if (isEMCal) {
ee1083d5 70 if (t.Contains("Detail")) {
71
72 cout<<"Rereading AliAnalysisEtMonteCarlo configuration file..."<<endl;
ee1083d5 73 fMCAnalysis = (AliAnalysisEmEtMonteCarlo *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(true,true)");
1815c8d7 74
ee1083d5 75 cout << "Instantiating AliAnalysisEmEtMonteCarlo class..."<< endl;
76 }
77 else if (fMCConfigFile.Length()) {
78 cout<<"Rereading AliAnalysisEtMonteCarloEmcal configuration file..."<<endl;
ee1083d5 79 fMCAnalysis = (AliAnalysisEtMonteCarloEmcal *) gInterpreter->ProcessLine("ConfigEtMonteCarlo()");
80 }
1815c8d7 81
ee1083d5 82 if (t.Contains("Detail")) {
83 cout<<"Rereading AliAnalysisEmEtReconstructed configuration file..."<<endl;
ee1083d5 84 fRecAnalysis = (AliAnalysisEmEtReconstructed *) gInterpreter->ProcessLine("ConfigEtReconstructed(true,true)");
85
86 }
87 else if (fRecoConfigFile.Length()) {
88 cout<<"Rereading AliAnalysisEtReconstructedEmcal configuration file..."<<endl;
ee1083d5 89 fRecAnalysis = (AliAnalysisEtReconstructedEmcal *) gInterpreter->ProcessLine("ConfigEtReconstructed()");
2fbf38ac 90 }
ee1083d5 91 }
92 else {
93 if (fMCConfigFile.Length()) {
94 cout<<"Rereading AliAnalysisEtMonteCarloPhos configuration file..."<<endl;
1815c8d7 95
ee1083d5 96 fMCAnalysis = (AliAnalysisEtMonteCarloPhos *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(false)");
97 cout << fMCAnalysis << endl;
98 }
1815c8d7 99
ee1083d5 100 if (fRecoConfigFile.Length()) {
101 cout<<"Rereading AliAnalysisEtReconstructedPhos configuration file..."<<endl;
ee1083d5 102 fRecAnalysis = (AliAnalysisEtReconstructedPhos *) gInterpreter->ProcessLine("ConfigEtReconstructed(false)");
2fbf38ac 103 }
ee1083d5 104 }
105 // Define input and output slots here
106 // Input slot #0 works with a TChain
107 DefineInput(0, TChain::Class());
108 // Output slot #1 writes into a TH1 container
1815c8d7 109
ee1083d5 110 DefineOutput(1, TList::Class());
1815c8d7 111
2fbf38ac 112}
d3d7bfe9 113AliAnalysisTaskTotEt::~AliAnalysisTaskTotEt() {//Destructor
ee1083d5 114 // fOutputList->Clear();
115 delete fRecAnalysis;
116 delete fMCAnalysis;
9a365626 117 delete fPIDResponse;
01b73fb0 118 //delete fSparseHistRecVsMc;
119 //delete fSparseRecVsMc;
951efd81 120}
2fbf38ac 121
122//________________________________________________________________________
123void AliAnalysisTaskTotEt::UserCreateOutputObjects()
124{
9a365626 125 //input hander
126 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
127 AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
128 if (!inputHandler) AliFatal("Input handler needed");
129
130 //pid response object
131 fPIDResponse=inputHandler->GetPIDResponse();
132 if (!fPIDResponse) AliError("PIDResponse object was not created");
133
ee1083d5 134 // Create histograms
135 // Called once
136 if (fMCAnalysis)
137 fMCAnalysis->CreateHistograms();
138 fRecAnalysis->CreateHistograms();
139 fOutputList = new TList;
140 fOutputList->SetOwner();
141 fRecAnalysis->FillOutputList(fOutputList);
142 if (fMCAnalysis)
143 fMCAnalysis->FillOutputList(fOutputList);
144 fHistEtRecvsEtMC = new TH2F("fHistEtRecvsEtMC", "Reconstructed E_{T} vs MC E_{T}", 1000, 0.000, 100, 1000, 0.0001, 100);
145 fHistEtRecOverEtMC = new TH2F("fHistEtRecOverEtMC", "Reconstructed E_{T} over MC E_{T} vs centrality", 1000, 0.00, 2.0, 11, -0.5, 10.5);
146 fHistDiffEtRecEtMCOverEtMC = new TH2F("fHistDiffEtRecEtMCOverEtMC", "fHistDiffEtRecEtMCOverEtMC", 10000, 0.0, 1000, 1000, -5, 5);
147 fOutputList->Add(fHistEtRecvsEtMC);
148 fOutputList->Add(fHistEtRecOverEtMC);
149 fOutputList->Add(fHistDiffEtRecEtMCOverEtMC);
150
151 Bool_t selectPrimaries=kTRUE;
152 if(fRecAnalysis->DataSet()==2010 || fRecAnalysis->DataSet()==20111||fRecAnalysis->DataSet()==2009){
153 if(fRecAnalysis->DataSet()==2010)cout<<"Setting track cuts for the 2010 p+p collisions at 7 TeV"<<endl;
154 else{
155 if(fRecAnalysis->DataSet()==2009){cout<<"Setting track cuts for the 2010 p+p collisions at 900 GeV"<<endl;}
156 else{cout<<"Setting track cuts for the 2011 p+p collisions at 2.76 TeV"<<endl;}
d3d7bfe9 157 }
ee1083d5 158 //cout<<"Warning: Have not set 2010 track cuts yet!!"<<endl;
159 fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
160 fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
161 fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
162 fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
163 //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
164 fEsdtrackCutsITS = AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
165 fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
166 }
167 if(fRecAnalysis->DataSet()==20100){
168 cout<<"Setting track cuts for the 2010 Pb+Pb collisions at 2.76 TeV"<<endl;
169 //cout<<"Warning: Have not set 2010 track cuts yet!!"<<endl;
170 fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
171 fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
172 fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
173 fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
174 //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
175 fEsdtrackCutsITS = AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
176 // fEsdtrackCutsITS = AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
177 fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
178 }
1815c8d7 179
ee1083d5 180
181 fOutputList->Add(fEsdtrackCutsITSTPC);
182 fOutputList->Add(fEsdtrackCutsTPC);
183 fOutputList->Add(fEsdtrackCutsITS);
184 if (fEsdtrackCutsITSTPC && fEsdtrackCutsTPC) {
185 fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
186 if (fMCAnalysis)
187 fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
188 fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
189 if (fMCAnalysis)
190 fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
191 fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
192 if (fMCAnalysis)
193 fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
194 //add ITS stuff!
195 }
196 else {
197 Printf("Error: no track cuts!");
198 }
1cc04082 199 PostData(1, fOutputList);
01b73fb0 200
2fbf38ac 201}
202
203//________________________________________________________________________
204void AliAnalysisTaskTotEt::UserExec(Option_t *)
cf6522d1 205{ // execute method
1815c8d7 206
ee1083d5 207 fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
208 if (!fESDEvent)
d3d7bfe9 209 {
ee1083d5 210 Printf("ERROR: Could not retrieve event");
211 return;
2fbf38ac 212 }
1815c8d7 213
ee1083d5 214 //Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber());
1815c8d7 215
ee1083d5 216 AliCentrality *cent = GetCentralityObject();
1815c8d7 217
ee1083d5 218 //if (res == 0 && cent)
219 //{
e50ad0d0 220 fRecAnalysis->SetCentralityObject(cent);
221 fRecAnalysis->AnalyseEvent(fESDEvent);
1815c8d7 222
e50ad0d0 223 AliMCEvent* mcEvent = MCEvent();
224 if (mcEvent)
225 {
226 fMCAnalysis->SetCentralityObject(cent);
227 fMCAnalysis->AnalyseEvent(mcEvent, fESDEvent);
228 //fMCAnalysis->AnalyseEvent(mcEvent);
229 }
230 if(fMCAnalysis)
231 {
f61cec2f 232 fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEt(), fMCAnalysis->GetTotNeutralEt());
233 if(fMCAnalysis->GetTotNeutralEt()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEt(), cent->GetCentralityClass10("V0M"));
234 if(fMCAnalysis->GetTotNeutralEt()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
ee1083d5 235 }
236 //}
237 // Post output data.
238 PostData(1, fOutputList);
2fbf38ac 239}
240
241//________________________________________________________________________
242void AliAnalysisTaskTotEt::Terminate(Option_t *)
243{
ee1083d5 244 // Draw result to the screen
245 // Called once at the end of the query
1815c8d7 246
ee1083d5 247 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
248 if (!fOutputList) {
249 printf("ERROR: Output list not available\n");
250 return;
251 }
2fbf38ac 252}
253
254
255