PWG4/totEt/AliAnalysisEt.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies
   3 //  Base class for ESD & MC analysis
   4 //  - reconstruction and MonteCarlo output
   5 // implementation file
   6 //
   7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
   8 //_________________________________________________________________________
   9
  10 #include "AliAnalysisEt.h"
  11 #include "TMath.h"
  12 #include "TList.h"
  13 #include "TH1F.h"
  14 #include "TH2F.h"
  15 #include <iostream>
  16 #include "AliAnalysisEtCuts.h"
  17 #include "AliVEvent.h"
  18 #include "TDatabasePDG.h"
  19 #include "Rtypes.h"
  20
  21 using namespace std;
  22 ClassImp(AliAnalysisEt);
  23
  24
  25 AliAnalysisEt::AliAnalysisEt() :
  26         fHistogramNameSuffix("")
  27         ,fCuts(0)
  28         ,fPdgDB(0)
  29         ,fPiPlusCode(0)
  30         ,fPiMinusCode(0)
  31         ,fKPlusCode(0)
  32         ,fKMinusCode(0)
  33         ,fProtonCode(0)
  34         ,fAntiProtonCode(0)
  35         ,fLambdaCode(0)
  36         ,fAntiLambdaCode(0)
  37         ,fK0SCode(0)
  38         ,fOmegaCode(0)
  39         ,fAntiOmegaCode(0)
  40         ,fXi0Code(0)
  41         ,fAntiXi0Code(0)
  42         ,fXiCode(0)
  43         ,fAntiXiCode(0)
  44         ,fSigmaCode(0)
  45         ,fAntiSigmaCode(0)
  46         ,fK0LCode(0)
  47         ,fNeutronCode(0)
  48         ,fAntiNeutronCode(0)
  49         ,fEPlusCode(0)
  50         ,fEMinusCode(0)
  51         ,fPionMass(0)
  52         ,fSumEt(0)
  53         ,fSumEtAcc(0)
  54         ,fTotEt(0)
  55         ,fTotEtAcc(0)
  56         ,fTotNeutralEt(0)
  57         ,fTotNeutralEtAcc(0)
  58         ,fTotChargedEt(0)
  59         ,fTotChargedEtAcc(0)
  60         ,fMultiplicity(0)
  61         ,fChargedMultiplicity(0)
  62         ,fNeutralMultiplicity(0)
  63         ,fEtaCutAcc(0)
  64         ,fPhiCutAccMin(0)
  65         ,fPhiCutAccMax(0)
  66         ,fDetectorRadius(0)
  67         ,fClusterEnergyCut(0)
  68         ,fSingleCellEnergyCut(0)
  69         ,fHistEt(0)
  70         ,fHistChargedEt(0)
  71         ,fHistNeutralEt(0)
  72         ,fHistEtAcc(0)
  73         ,fHistChargedEtAcc(0)
  74         ,fHistNeutralEtAcc(0)
  75         ,fHistMult(0)
  76         ,fHistChargedMult(0)
  77         ,fHistNeutralMult(0)
  78         ,fHistPhivsPtPos(0)
  79         ,fHistPhivsPtNeg(0)
  80         ,fHistBaryonEt(0)
  81         ,fHistAntiBaryonEt(0)
  82         ,fHistMesonEt(0)
  83         ,fHistBaryonEtAcc(0)
  84         ,fHistAntiBaryonEtAcc(0)
  85         ,fHistMesonEtAcc(0)
  86         ,fHistTMDeltaR(0)
  87 {
  88 }
  89
  90 AliAnalysisEt::~AliAnalysisEt()
  91 {
  92
  93 }
  94
  95 void AliAnalysisEt::FillOutputList(TList *list)
  96 { // histograms to be added to output
  97     list->Add(fHistEt);
  98     list->Add(fHistChargedEt);
  99     list->Add(fHistNeutralEt);
 100
 101     list->Add(fHistEtAcc);
 102     list->Add(fHistChargedEtAcc);
 103     list->Add(fHistNeutralEtAcc);
 104
 105     list->Add(fHistMult);
 106     list->Add(fHistChargedMult);
 107     list->Add(fHistNeutralMult);
 108
 109     list->Add(fHistPhivsPtPos);
 110     list->Add(fHistPhivsPtNeg);
 111
 112     list->Add(fHistBaryonEt);
 113     list->Add(fHistAntiBaryonEt);
 114     list->Add(fHistMesonEt);
 115
 116     list->Add(fHistBaryonEtAcc);
 117     list->Add(fHistAntiBaryonEtAcc);
 118     list->Add(fHistMesonEtAcc);
 119
 120     list->Add(fHistTMDeltaR);
 121 }
 122
 123 void AliAnalysisEt::Init()
 124 {// clear variables, set up cuts and PDG info
 125   ResetEventValues();
 126 }
 127
 128 void AliAnalysisEt::CreateHistograms()
 129 { // create histograms..
 130   // histogram binning for E_T, p_T and Multiplicity: defaults for p+p
 131   Int_t nbinsEt = 1000;
 132   Double_t minEt = 0.0001;
 133   Double_t maxEt = 100;
 134   Int_t nbinsPt = 200;
 135   Double_t minPt = 0;
 136   Double_t maxPt = 20;
 137   Int_t nbinsMult = 200;
 138   Double_t minMult = -0.5; // offset -0.5 to have integer bins centered around 0
 139   Double_t maxMult = nbinsMult + minMult; // 1 bin per integer value
 140
 141     TString histname = "fHistEt" + fHistogramNameSuffix;
 142     fHistEt = new TH1F(histname.Data(), "Total E_{T} Distribution", nbinsEt, minEt, maxEt);
 143     fHistEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 144     fHistEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 145
 146     histname = "fHistChargedEt" + fHistogramNameSuffix;
 147     fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", nbinsEt, minEt, maxEt);
 148     fHistChargedEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 149     fHistChargedEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 150
 151     histname = "fHistNeutralEt" + fHistogramNameSuffix;
 152     fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", nbinsEt, minEt, maxEt);
 153     fHistNeutralEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 154     fHistNeutralEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 155
 156     histname = "fHistEtAcc" + fHistogramNameSuffix;
 157     fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
 158     fHistEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 159     fHistEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 160
 161     histname = "fHistChargedEtAcc" + fHistogramNameSuffix;
 162     fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
 163     fHistChargedEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 164     fHistChargedEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 165
 166     histname = "fHistNeutralEtAcc" + fHistogramNameSuffix;
 167     fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
 168     fHistNeutralEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
 169     fHistNeutralEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 170     std::cout << histname << std::endl;
 171     histname = "fHistMult" + fHistogramNameSuffix;
 172     fHistMult = new TH1F(histname.Data(), "Total Multiplicity", nbinsMult, minMult, maxMult);
 173     fHistMult->GetXaxis()->SetTitle("N");
 174     fHistMult->GetYaxis()->SetTitle("Multiplicity");
 175
 176     histname = "fHistChargedMult" + fHistogramNameSuffix;
 177     fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", nbinsMult, minMult, maxMult);
 178     fHistChargedMult->GetXaxis()->SetTitle("N");
 179     fHistChargedMult->GetYaxis()->SetTitle("Multiplicity");
 180
 181     histname = "fHistNeutralMult" + fHistogramNameSuffix;
 182     fHistNeutralMult = new TH1F(histname.Data(), "Neutral Multiplicity", nbinsMult, minMult, maxMult);
 183     fHistNeutralMult->GetXaxis()->SetTitle("N");
 184     fHistNeutralMult->GetYaxis()->SetTitle("Multiplicity");
 185
 186     histname = "fHistPhivsPtPos" + fHistogramNameSuffix;
 187     fHistPhivsPtPos = new TH2F(histname.Data(), "Phi vs pT of positively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
 188
 189     histname = "fHistPhivsPtNeg" + fHistogramNameSuffix;
 190     fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
 191
 192     histname = "fHistBaryonEt" + fHistogramNameSuffix;
 193     fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons",  nbinsEt, minEt, maxEt);
 194
 195     histname = "fHistAntiBaryonEt" + fHistogramNameSuffix;
 196     fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons",  nbinsEt, minEt, maxEt);
 197
 198     histname = "fHistMesonEt" + fHistogramNameSuffix;
 199     fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons",  nbinsEt, minEt, maxEt);
 200
 201     histname = "fHistBaryonEtAcc" + fHistogramNameSuffix;
 202     fHistBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
 203
 204     histname = "fHistAntiBaryonEtAcc" + fHistogramNameSuffix;
 205     fHistAntiBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for anti baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
 206
 207     histname = "fHistMesonEtAcc" + fHistogramNameSuffix;
 208     fHistMesonEtAcc = new TH1F(histname.Data(), "E_{T} for mesons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
 209
 210     //
 211     histname = "fHistTMDeltaR" + fHistogramNameSuffix;
 212     fHistTMDeltaR = new TH1F(histname.Data(), "#Delta R for calorimeter clusters", 200, 0, 50);
 213
 214 }
 215
 216 void AliAnalysisEt::FillHistograms()
 217 { // fill histograms..
 218     fHistEt->Fill(fTotEt);
 219     fHistChargedEt->Fill(fTotChargedEt);
 220     fHistNeutralEt->Fill(fTotNeutralEt);
 221
 222     fHistEtAcc->Fill(fTotEtAcc);
 223     fHistChargedEtAcc->Fill(fTotChargedEtAcc);
 224     fHistNeutralEtAcc->Fill(fTotNeutralEtAcc);
 225
 226     fHistMult->Fill(fMultiplicity);
 227     fHistChargedMult->Fill(fChargedMultiplicity);
 228     fHistNeutralMult->Fill(fNeutralMultiplicity);
 229
 230     /* // DS commented out non-fills to prevent compilation warnings
 231     fHistPhivsPtPos;
 232     fHistPhivsPtNeg;
 233
 234     fHistBaryonEt;
 235     fHistAntiBaryonEt;
 236     fHistMesonEt;
 237
 238     fHistBaryonEtAcc;
 239     fHistAntiBaryonEtAcc;
 240     fHistMesonEtAcc;
 241
 242     fHistTMDeltaR;
 243     */
 244 }
 245
 246 Int_t AliAnalysisEt::AnalyseEvent(AliVEvent *event)
 247 { //this line is basically here to eliminate a compiler warning that event is not used.  Making it a virtual function did not work with the plugin.
 248   cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
 249   ResetEventValues();
 250   return 0;
 251 }
 252
 253 void AliAnalysisEt::ResetEventValues()
 254 { // clear
 255   fTotEt = 0;
 256   fTotEtAcc = 0;
 257   fTotNeutralEt = 0;
 258   fTotNeutralEtAcc = 0;
 259   fTotChargedEt  = 0;
 260   fTotChargedEtAcc = 0;
 261   fMultiplicity = 0;
 262   fChargedMultiplicity = 0;
 263   fNeutralMultiplicity = 0;
 264
 265   if (!fCuts || !fPdgDB || fPiPlusCode==0) { // some Init's needed
 266     cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
 267     if (!fCuts) {
 268       cout << " setting up Cuts " << endl;
 269       fCuts = new AliAnalysisEtCuts();
 270     }
 271     if(!fPdgDB) {
 272       cout << " setting up PdgDB " << endl;
 273       fPdgDB = new TDatabasePDG();
 274     }
 275
 276     if (fPiPlusCode==0) {
 277       SetParticleCodes();
 278     }
 279   }
 280   return;
 281 }
 282
 283
 284 void AliAnalysisEt::SetParticleCodes()
 285 { // set PDG info
 286   fPionMass = fPdgDB->GetParticle("pi+")->Mass();
 287   fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
 288   fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
 289   fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
 290   fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
 291   fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
 292   fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
 293   fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
 294   fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
 295   fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
 296   fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
 297   fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
 298   fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
 299   fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
 300   fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
 301   fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
 302   fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
 303   fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
 304   fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
 305   fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
 306   fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
 307   fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
 308   fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
 309   cout << "Resetting Codes: Pion " << fPiPlusCode
 310        << "," << fPiMinusCode
 311        << " Kaon " << fKPlusCode
 312        << "," << fKMinusCode << endl;
 313 }
 314