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