1 //_________________________________________________________________________
2 // Utility Class for transverse energy studies
4 // - reconstruction and MC output
7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8 //_________________________________________________________________________
9 //Necessary to read config macros
12 #include <TInterpreter.h>
18 #include "THnSparse.h"
20 #include "AliESDEvent.h"
21 #include "AliMCEvent.h"
22 #include "AliESDtrackCuts.h"
24 #include "AliAnalysisTaskTotEt.h"
25 #include "AliAnalysisEtReconstructedPhos.h"
26 #include "AliAnalysisEtReconstructedEmcal.h"
27 #include "AliAnalysisEtMonteCarloPhos.h"
28 #include "AliAnalysisEtMonteCarloEmcal.h"
29 #include "AliAnalysisEmEtMonteCarlo.h"
30 #include "AliAnalysisEmEtReconstructed.h"
31 #include "AliAnalysisManager.h"
32 #include "AliPIDResponse.h"
33 #include "AliTPCPIDResponse.h"
34 #include "AliInputEventHandler.h"
37 #include <AliCentrality.h>
41 ClassImp(AliAnalysisTaskTotEt)
43 //________________________________________________________________________
44 AliAnalysisTaskTotEt::AliAnalysisTaskTotEt(const char *name, Bool_t isMc) :
45 AliAnalysisTaskTransverseEnergy(name, isMc)
49 // ,fSparseHistRecVsMc(0)
53 // select if we should use EMCal or PHOS class
54 // PHOS by default, EMCal if name string contains EMC
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()");
67 //cout<<__FILE__<<" My name is "<<name<<endl;
70 if (t.Contains("Detail")) {
72 cout<<"Rereading AliAnalysisEtMonteCarlo configuration file..."<<endl;
73 fMCAnalysis = (AliAnalysisEmEtMonteCarlo *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(true,true)");
75 cout << "Instantiating AliAnalysisEmEtMonteCarlo class..."<< endl;
77 else if (fMCConfigFile.Length()) {
78 cout<<"Rereading AliAnalysisEtMonteCarloEmcal configuration file..."<<endl;
79 fMCAnalysis = (AliAnalysisEtMonteCarloEmcal *) gInterpreter->ProcessLine("ConfigEtMonteCarlo()");
82 if (t.Contains("Detail")) {
83 cout<<"Rereading AliAnalysisEmEtReconstructed configuration file..."<<endl;
84 fRecAnalysis = (AliAnalysisEmEtReconstructed *) gInterpreter->ProcessLine("ConfigEtReconstructed(true,true)");
87 else if (fRecoConfigFile.Length()) {
88 cout<<"Rereading AliAnalysisEtReconstructedEmcal configuration file..."<<endl;
89 fRecAnalysis = (AliAnalysisEtReconstructedEmcal *) gInterpreter->ProcessLine("ConfigEtReconstructed()");
93 if (fMCConfigFile.Length()) {
94 cout<<"Rereading AliAnalysisEtMonteCarloPhos configuration file..."<<endl;
96 fMCAnalysis = (AliAnalysisEtMonteCarloPhos *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(false)");
97 cout << fMCAnalysis << endl;
100 if (fRecoConfigFile.Length()) {
101 cout<<"Rereading AliAnalysisEtReconstructedPhos configuration file..."<<endl;
102 fRecAnalysis = (AliAnalysisEtReconstructedPhos *) gInterpreter->ProcessLine("ConfigEtReconstructed(false)");
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
110 DefineOutput(1, TList::Class());
113 AliAnalysisTaskTotEt::~AliAnalysisTaskTotEt() {//Destructor
114 // fOutputList->Clear();
118 //delete fSparseHistRecVsMc;
119 //delete fSparseRecVsMc;
122 //________________________________________________________________________
123 void AliAnalysisTaskTotEt::UserCreateOutputObjects()
126 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
127 AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
128 if (!inputHandler) AliFatal("Input handler needed");
130 //pid response object
131 fPIDResponse=inputHandler->GetPIDResponse();
132 if (!fPIDResponse) AliError("PIDResponse object was not created");
137 fMCAnalysis->CreateHistograms();
138 fRecAnalysis->CreateHistograms();
139 fOutputList = new TList;
140 fOutputList->SetOwner();
141 fRecAnalysis->FillOutputList(fOutputList);
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);
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;
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;}
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");
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");
181 fOutputList->Add(fEsdtrackCutsITSTPC);
182 fOutputList->Add(fEsdtrackCutsTPC);
183 fOutputList->Add(fEsdtrackCutsITS);
184 if (fEsdtrackCutsITSTPC && fEsdtrackCutsTPC) {
185 fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
187 fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
188 fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
190 fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
191 fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
193 fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
197 Printf("Error: no track cuts!");
199 PostData(1, fOutputList);
203 //________________________________________________________________________
204 void AliAnalysisTaskTotEt::UserExec(Option_t *)
207 fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
210 Printf("ERROR: Could not retrieve event");
214 //Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber());
216 AliCentrality *cent = GetCentralityObject();
218 //if (res == 0 && cent)
220 fRecAnalysis->SetCentralityObject(cent);
221 fRecAnalysis->AnalyseEvent(fESDEvent);
223 AliMCEvent* mcEvent = MCEvent();
226 fMCAnalysis->SetCentralityObject(cent);
227 fMCAnalysis->AnalyseEvent(mcEvent, fESDEvent);
228 //fMCAnalysis->AnalyseEvent(mcEvent);
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());
238 PostData(1, fOutputList);
241 //________________________________________________________________________
242 void AliAnalysisTaskTotEt::Terminate(Option_t *)
244 // Draw result to the screen
245 // Called once at the end of the query
247 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
249 printf("ERROR: Output list not available\n");