]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/totEt/AliAnalysisTaskTotEt.cxx
add default options for kepler gpu family
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisTaskTotEt.cxx
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                //_________________________________________________________________________
9                //Necessary to read config macros
10 #include <TROOT.h>
11 #include <TSystem.h>
12 #include <TInterpreter.h>
13
14 #include "TChain.h"
15 #include "TList.h"
16 #include "TFile.h"
17 #include "TH2F.h"
18 #include "THnSparse.h"
19
20 #include "AliESDEvent.h"
21 #include "AliMCEvent.h"
22 #include "AliESDtrackCuts.h"
23
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"
35
36 #include <iostream>
37 #include <AliCentrality.h>
38
39   using namespace std;
40
41 ClassImp(AliAnalysisTaskTotEt)
42
43 //________________________________________________________________________
44   AliAnalysisTaskTotEt::AliAnalysisTaskTotEt(const char *name, Bool_t isMc) :
45     AliAnalysisTaskTransverseEnergy(name, isMc)
46     ,fPIDResponse(0)
47     ,fRecAnalysis(0)
48     ,fMCAnalysis(0)
49                                             // ,fSparseHistRecVsMc(0)
50                                             //,fSparseRecVsMc(0)
51 {
52   // Constructor
53   // select if we should use EMCal or PHOS class
54   // PHOS by default, EMCal if name string contains EMC
55   TString t(name);
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;
68   //t.ToUpper();
69   if (isEMCal) {
70     if (t.Contains("Detail")) {
71
72       cout<<"Rereading AliAnalysisEtMonteCarlo configuration file..."<<endl;
73       fMCAnalysis = (AliAnalysisEmEtMonteCarlo *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(true,true)");
74                         
75       cout << "Instantiating AliAnalysisEmEtMonteCarlo class..."<< endl;
76     }
77     else if (fMCConfigFile.Length()) {
78       cout<<"Rereading AliAnalysisEtMonteCarloEmcal configuration file..."<<endl;
79       fMCAnalysis = (AliAnalysisEtMonteCarloEmcal *) gInterpreter->ProcessLine("ConfigEtMonteCarlo()");
80     }
81                 
82     if (t.Contains("Detail")) {
83       cout<<"Rereading AliAnalysisEmEtReconstructed configuration file..."<<endl;
84       fRecAnalysis = (AliAnalysisEmEtReconstructed *) gInterpreter->ProcessLine("ConfigEtReconstructed(true,true)");
85
86     }
87     else if (fRecoConfigFile.Length()) {
88       cout<<"Rereading AliAnalysisEtReconstructedEmcal configuration file..."<<endl;
89       fRecAnalysis = (AliAnalysisEtReconstructedEmcal *) gInterpreter->ProcessLine("ConfigEtReconstructed()");
90     }
91   }
92   else {
93     if (fMCConfigFile.Length()) {
94       cout<<"Rereading AliAnalysisEtMonteCarloPhos configuration file..."<<endl;
95                         
96       fMCAnalysis = (AliAnalysisEtMonteCarloPhos *) gInterpreter->ProcessLine("ConfigEtMonteCarlo(false)");
97       cout << fMCAnalysis << endl;
98     }
99                 
100     if (fRecoConfigFile.Length()) {
101       cout<<"Rereading AliAnalysisEtReconstructedPhos configuration file..."<<endl;
102       fRecAnalysis = (AliAnalysisEtReconstructedPhos *) gInterpreter->ProcessLine("ConfigEtReconstructed(false)");
103     }
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
109         
110   DefineOutput(1, TList::Class());
111         
112 }
113 AliAnalysisTaskTotEt::~AliAnalysisTaskTotEt() {//Destructor
114   //    fOutputList->Clear();
115   delete fRecAnalysis;
116   delete fMCAnalysis;
117   delete fPIDResponse;
118   //delete fSparseHistRecVsMc;
119   //delete fSparseRecVsMc;
120 }
121
122 //________________________________________________________________________
123 void AliAnalysisTaskTotEt::UserCreateOutputObjects()
124 {
125   //input hander
126   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
127   if (!man) {
128     AliFatal("Analysis manager needed");
129     return;
130   }
131
132   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
133   if (!inputHandler) {
134     AliFatal("Input handler needed");
135     return;
136   }
137
138   //pid response object
139   fPIDResponse=inputHandler->GetPIDResponse();
140   if (!fPIDResponse) AliError("PIDResponse object was not created");
141
142   // Create histograms
143   // Called once
144   if (fMCAnalysis)
145     fMCAnalysis->CreateHistograms();
146   fRecAnalysis->CreateHistograms();
147   fOutputList = new TList;
148   fOutputList->SetOwner();
149   fRecAnalysis->FillOutputList(fOutputList);
150   if (fMCAnalysis)
151     fMCAnalysis->FillOutputList(fOutputList);
152   fHistEtRecvsEtMC = new TH2F("fHistEtRecvsEtMC", "Reconstructed E_{T} vs MC E_{T}", 1000, 0.000, 100, 1000, 0.0001, 100);
153   fHistEtRecOverEtMC = new TH2F("fHistEtRecOverEtMC", "Reconstructed E_{T} over MC E_{T} vs centrality", 1000, 0.00, 2.0, 11, -0.5, 10.5); 
154   fHistDiffEtRecEtMCOverEtMC = new TH2F("fHistDiffEtRecEtMCOverEtMC", "fHistDiffEtRecEtMCOverEtMC", 10000, 0.0, 1000, 1000, -5, 5); 
155   fOutputList->Add(fHistEtRecvsEtMC);
156   fOutputList->Add(fHistEtRecOverEtMC);
157   fOutputList->Add(fHistDiffEtRecEtMCOverEtMC);
158
159   Bool_t selectPrimaries=kTRUE;
160   if(fRecAnalysis->DataSet()==2010 || fRecAnalysis->DataSet()==20111||fRecAnalysis->DataSet()==2009 || fRecAnalysis->DataSet()==2012 || fRecAnalysis->DataSet()==2013){
161     if(fRecAnalysis->DataSet()==2010)cout<<"Setting track cuts for the 2010 p+p collisions at 7 TeV"<<endl;
162     else{
163       if(fRecAnalysis->DataSet()==2012)cout<<"Setting track cuts for the 2012 p+p collisions at 8 TeV"<<endl;
164       else{
165         if(fRecAnalysis->DataSet()==2013)cout<<"Setting track cuts for the 2013 p+Pb collisions at 5 TeV"<<endl;
166         else{
167           if(fRecAnalysis->DataSet()==2009){cout<<"Setting track cuts for the 2010 p+p collisions at 900 GeV"<<endl;}
168           else{cout<<"Setting track cuts for the 2011 p+p collisions at 2.76 TeV"<<endl;}
169         }
170       }
171     }
172     //cout<<"Warning:  Have not set 2010 track cuts yet!!"<<endl;
173     fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
174     fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
175     fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
176     fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
177     //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
178     fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
179     fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
180   }
181   if(fRecAnalysis->DataSet()==20100){
182     cout<<"Setting track cuts for the 2010 Pb+Pb collisions at 2.76 TeV"<<endl;
183     //cout<<"Warning:  Have not set 2010 track cuts yet!!"<<endl;
184     fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
185     fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
186     fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
187     fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
188     //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
189     fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
190     // fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
191    fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
192   }
193         
194         
195   fOutputList->Add(fEsdtrackCutsITSTPC);
196   fOutputList->Add(fEsdtrackCutsTPC);
197   fOutputList->Add(fEsdtrackCutsITS);
198   if (fEsdtrackCutsITSTPC && fEsdtrackCutsTPC) {
199     fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
200     if (fMCAnalysis)
201       fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
202     fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
203     if (fMCAnalysis)
204       fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
205     fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
206     if (fMCAnalysis)
207       fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
208     //add ITS stuff!
209   }
210   else {
211     Printf("Error: no track cuts!");
212   }
213  PostData(1, fOutputList);
214
215 }
216
217 //________________________________________________________________________
218 void AliAnalysisTaskTotEt::UserExec(Option_t *)
219 { // execute method
220         
221   fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
222   if (!fESDEvent)
223     {
224       Printf("ERROR: Could not retrieve event");
225       return;
226     }
227         
228   //Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber());
229         
230   AliCentrality *cent = GetCentralityObject();
231         
232   //if (res == 0 && cent)
233   //{
234   fRecAnalysis->SetCentralityObject(cent);
235   fRecAnalysis->AnalyseEvent(fESDEvent);
236                 
237   AliMCEvent* mcEvent = MCEvent();
238   if (mcEvent)
239     {
240       fMCAnalysis->SetCentralityObject(cent);
241       fMCAnalysis->AnalyseEvent(mcEvent, fESDEvent);
242       //fMCAnalysis->AnalyseEvent(mcEvent);
243     }
244   if(fMCAnalysis)
245     {
246       fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEt(), fMCAnalysis->GetTotNeutralEt());
247       if(fMCAnalysis->GetTotNeutralEt()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEt(), cent->GetCentralityClass10("V0M"));
248       if(fMCAnalysis->GetTotNeutralEt()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
249     }
250   //}
251   // Post output data.
252   PostData(1, fOutputList);
253 }
254
255 //________________________________________________________________________
256 void AliAnalysisTaskTotEt::Terminate(Option_t *)
257 {
258   // Draw result to the screen
259   // Called once at the end of the query
260         
261   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
262   if (!fOutputList) {
263     printf("ERROR: Output list not available\n");
264     return;
265   }
266 }
267
268
269