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