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