PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassScale.cxx

   1 //
   2 // rho mass scale task
   3 // task to estimate scale factor for rho_m
   4 //
   5 // Author: M.Verweij
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <THnSparse.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TProfile.h>
  15 #include <TChain.h>
  16 #include <TSystem.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19
  20 #include "AliAnalysisManager.h"
  21 // #include "AliAODMCHeader.h"
  22 // #include "AliMCEvent.h"
  23 // #include "AliGenPythiaEventHeader.h"
  24 // #include "AliAODEvent.h"
  25 #include "AliLog.h"
  26 #include "AliVCluster.h"
  27 #include "AliVTrack.h"
  28 #include "AliEmcalJet.h"
  29 #include "AliRhoParameter.h"
  30 #include "AliEmcalParticle.h"
  31 #include "AliJetContainer.h"
  32 #include "AliParticleContainer.h"
  33
  34 #include "AliAnalysisTaskRhoMassScale.h"
  35
  36 ClassImp(AliAnalysisTaskRhoMassScale)
  37
  38 //________________________________________________________________________
  39 AliAnalysisTaskRhoMassScale::AliAnalysisTaskRhoMassScale() :
  40   AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoMassScale", kTRUE),
  41   fContainerNeutral(0),
  42   fContainerCharged(1),
  43   fRhoMNeutralName(""),
  44   fRhoMChargedEmcalName(""),
  45   fRhoMCharged2xEmcalName(""),
  46   fRhoMNeutral(0),
  47   fRhoMChargedEmcal(0),
  48   fRhoMCharged2xEmcal(0),
  49   fHistScaleEmcalvsCent(0),
  50   fHistScale2EmcalvsCent(0),
  51   fHistDeltaScale2EmcalvsCent(0),
  52   fHistScaleEmcalvsMult(0),
  53   fHistScale2EmcalvsMult(0),
  54   fHistDeltaScale2EmcalvsMult(0)
  55 {
  56   // Default constructor.
  57
  58   SetMakeGeneralHistograms(kTRUE);
  59 }
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskRhoMassScale::AliAnalysisTaskRhoMassScale(const char *name) :
  63   AliAnalysisTaskEmcalJet(name, kTRUE),
  64   fContainerNeutral(0),
  65   fContainerCharged(1),
  66   fRhoMNeutralName(""),
  67   fRhoMChargedEmcalName(""),
  68   fRhoMCharged2xEmcalName(""),
  69   fRhoMNeutral(0),
  70   fRhoMChargedEmcal(0),
  71   fRhoMCharged2xEmcal(0),
  72   fHistScaleEmcalvsCent(0),
  73   fHistScale2EmcalvsCent(0),
  74   fHistDeltaScale2EmcalvsCent(0),
  75   fHistScaleEmcalvsMult(0),
  76   fHistScale2EmcalvsMult(0),
  77   fHistDeltaScale2EmcalvsMult(0)
  78 {
  79   // Standard constructor.
  80
  81   SetMakeGeneralHistograms(kTRUE);
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskRhoMassScale::~AliAnalysisTaskRhoMassScale()
  86 {
  87   // Destructor.
  88 }
  89
  90 //________________________________________________________________________
  91 void AliAnalysisTaskRhoMassScale::UserCreateOutputObjects()
  92 {
  93   // Create user output.
  94
  95   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  96
  97   Bool_t oldStatus = TH1::AddDirectoryStatus();
  98   TH1::AddDirectory(kFALSE);
  99
 100   //Create histograms
 101   TString histName = "";
 102   TString histTitle = "";
 103
 104   histName = "fHistScaleEmcalvsCent";
 105   histTitle = TString::Format("%s;Centrality;s_{EMC}",histName.Data());
 106   fHistScaleEmcalvsCent= new TH2F(histName.Data(),histTitle.Data(), 101, -1, 100, 500, 0, 5);
 107   fOutput->Add(fHistScaleEmcalvsCent);
 108
 109   histName = "fHistScale2EmcalvsCent";
 110   histTitle = TString::Format("%s;Centrality;s_{2 #times EMC}",histName.Data());
 111   fHistScale2EmcalvsCent = new TH2F(histName.Data(),histTitle.Data(), 101, -1, 100, 500, 0, 5);
 112   fOutput->Add(fHistScale2EmcalvsCent);
 113
 114   histName = "fHistDeltaScale2EmcalvsCent";
 115   histTitle = TString::Format("%s;Centrality;s_{2 #times EMC}-s_{EMC}",histName.Data());
 116   fHistDeltaScale2EmcalvsCent = new TH2F(histName.Data(),histTitle.Data(), 101, -1, 100, 500, -2.5, 2.5);
 117   fOutput->Add(fHistDeltaScale2EmcalvsCent);
 118
 119   histName = "fHistScaleEmcalvsMult";
 120   histTitle = TString::Format("%s;#it{N}_{track};s_{EMC}",histName.Data());
 121   fHistScaleEmcalvsMult= new TH2F(histName.Data(),histTitle.Data(), 800, 0, 4000, 500, 0, 5);
 122   fOutput->Add(fHistScaleEmcalvsMult);
 123
 124   histName = "fHistScale2EmcalvsMult";
 125   histTitle = TString::Format("%s;#it{N}_{track};s_{2 #times EMC}",histName.Data());
 126   fHistScale2EmcalvsMult = new TH2F(histName.Data(),histTitle.Data(), 800, 0, 4000, 500, 0, 5);
 127   fOutput->Add(fHistScale2EmcalvsMult);
 128
 129   histName = "fHistDeltaScale2EmcalvsMult";
 130   histTitle = TString::Format("%s;#it{N}_{track};s_{2 #times EMC}-s_{EMC}",histName.Data());
 131   fHistDeltaScale2EmcalvsMult = new TH2F(histName.Data(),histTitle.Data(), 800, 0, 4000, 500, -2.5, 2.5);
 132   fOutput->Add(fHistDeltaScale2EmcalvsMult);
 133
 134   // =========== Switch on Sumw2 for all histos ===========
 135   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 136     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 137     if (h1){
 138       h1->Sumw2();
 139       continue;
 140     }
 141     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 142     if(hn)hn->Sumw2();
 143   }
 144
 145   TH1::AddDirectory(oldStatus);
 146
 147   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 148 }
 149
 150 //________________________________________________________________________
 151 Bool_t AliAnalysisTaskRhoMassScale::Run()
 152 {
 153   // Run analysis code here, if needed. It will be executed before FillHistograms().
 154   return kTRUE;
 155 }
 156
 157 //________________________________________________________________________
 158 Bool_t AliAnalysisTaskRhoMassScale::FillHistograms()
 159 {
 160   // Fill histograms.
 161
 162   Double_t rhomNe        = fRhoMNeutral->GetVal();
 163   Double_t rhomChEmcal   = fRhoMChargedEmcal->GetVal();
 164   Double_t rhomCh2xEmcal = fRhoMCharged2xEmcal->GetVal();
 165
 166   Double_t scale = -1.; Double_t scale2 = -1.;
 167   if(rhomChEmcal>0.)   scale = (rhomNe+rhomChEmcal)/rhomChEmcal;
 168   if(rhomCh2xEmcal>0.) scale2 = (rhomNe+rhomChEmcal)/rhomCh2xEmcal;
 169
 170   fHistScaleEmcalvsCent->Fill(fCent,scale);
 171   fHistScale2EmcalvsCent->Fill(fCent,scale2);
 172   fHistDeltaScale2EmcalvsCent->Fill(fCent,scale2-scale);
 173
 174   Int_t mult = -1;
 175   if(GetParticleContainer(0))
 176     mult = GetParticleContainer(0)->GetNAcceptedParticles();
 177
 178   fHistScaleEmcalvsMult->Fill(mult,scale);
 179   fHistScale2EmcalvsMult->Fill(mult,scale2);
 180   fHistDeltaScale2EmcalvsMult->Fill(mult,scale2-scale);
 181
 182   return kTRUE;
 183 }
 184
 185 //________________________________________________________________________
 186 Bool_t AliAnalysisTaskRhoMassScale::RetrieveEventObjects() {
 187   //
 188   // retrieve event objects
 189   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 190     return kFALSE;
 191
 192   if (!fRhoMNeutralName.IsNull() && !fRhoMNeutral) { // get rho_m from the event
 193     fRhoMNeutral = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoMNeutralName));
 194     if (!fRhoMNeutral) {
 195       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoMNeutralName.Data()));
 196       fInitialized = kFALSE;
 197       return kFALSE;
 198     }
 199   }
 200
 201   if (!fRhoMChargedEmcalName.IsNull() && !fRhoMChargedEmcal) { // get rho_m from the event
 202     fRhoMChargedEmcal = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoMChargedEmcalName));
 203     if (!fRhoMChargedEmcal) {
 204       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoMChargedEmcalName.Data()));
 205       fInitialized = kFALSE;
 206       return kFALSE;
 207     }
 208   }
 209
 210   if (!fRhoMCharged2xEmcalName.IsNull() && !fRhoMCharged2xEmcal) { // get rho_m from the event
 211     fRhoMCharged2xEmcal = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoMCharged2xEmcalName));
 212     if (!fRhoMCharged2xEmcal) {
 213       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoMCharged2xEmcalName.Data()));
 214       fInitialized = kFALSE;
 215       return kFALSE;
 216     }
 217   }
 218
 219   return kTRUE;
 220 }
 221
 222 //_______________________________________________________________________
 223 void AliAnalysisTaskRhoMassScale::Terminate(Option_t *)
 224 {
 225   // Called once at the end of the analysis.
 226 }
 227