PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaAnalysisMCSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TParticle.h>
   9 #include <TParticlePDG.h>
  10 #include <TVector3.h>
  11 #include <TH1F.h>
  12 #include <TH3F.h>
  13 #include <TTree.h>
  14
  15 #include <AliLog.h>
  16 #include <AliGenEventHeader.h>
  17 #include <AliHeader.h>
  18
  19 #include "dNdEtaAnalysis.h"
  20
  21
  22 ClassImp(AlidNdEtaAnalysisMCSelector)
  23
  24 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
  25   AlidNdEtaAnalysisSelector(),
  26   fVertex(0),
  27   fPartEta(0),
  28   fEvents(0)
  29 {
  30   //
  31   // Constructor. Initialization of pointers
  32   //
  33 }
  34
  35 AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
  36 {
  37   //
  38   // Destructor
  39   //
  40 }
  41
  42 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
  43 {
  44    AlidNdEtaAnalysisSelector::Init(tree);
  45
  46   tree->SetBranchStatus("ESD", 0);
  47
  48   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
  49   fPartEta = new TH1F("dndeta_check", "dndeta_check", 120, -6, 6);
  50   fPartEta->Sumw2();
  51 }
  52
  53 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
  54 {
  55   //
  56
  57   if (AliSelector::Process(entry) == kFALSE)
  58     return kFALSE;
  59
  60   TTree* particleTree = GetKinematics();
  61   if (!particleTree)
  62   {
  63     AliDebug(AliLog::kError, "Kinematics not available");
  64     return kFALSE;
  65   }
  66
  67   AliHeader* header = GetHeader();
  68   if (!header)
  69   {
  70     AliDebug(AliLog::kError, "Header not available");
  71     return kFALSE;
  72   }
  73
  74   // get the MC vertex
  75   AliGenEventHeader* genHeader = header->GenEventHeader();
  76
  77   TArrayF vtxMC(3);
  78   genHeader->PrimaryVertex(vtxMC);
  79
  80   particleTree->SetBranchStatus("*", 0);
  81   particleTree->SetBranchStatus("fDaughter[2]", 1);
  82   particleTree->SetBranchStatus("fPdgCode", 1);
  83   particleTree->SetBranchStatus("fPx", 1);
  84   particleTree->SetBranchStatus("fPy", 1);
  85   particleTree->SetBranchStatus("fPz", 1);
  86   particleTree->SetBranchStatus("fVx", 1);
  87   particleTree->SetBranchStatus("fVy", 1);
  88   particleTree->SetBranchStatus("fVz", 1);
  89
  90   TParticle* particle = 0;
  91   particleTree->SetBranchAddress("Particles", &particle);
  92
  93   Int_t nPrim  = header->GetNprimary();
  94   Int_t nTotal = header->GetNtrack();
  95
  96   for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
  97   {
  98     particleTree->GetEntry(i_mc);
  99
 100     if (!particle)
 101       continue;
 102
 103     if (IsPrimaryCharged(particle, nPrim) == kFALSE)
 104       continue;
 105
 106     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", i_mc, particle->GetUniqueID()));
 107
 108     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta());
 109     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 110
 111     fPartEta->Fill(particle->Eta());
 112   }
 113   fdNdEtaAnalysis->FillEvent(vtxMC[2]);
 114
 115   ++fEvents;
 116
 117   return kTRUE;
 118 }
 119
 120 void AlidNdEtaAnalysisMCSelector::Terminate()
 121 {
 122   AlidNdEtaAnalysisSelector::Terminate();
 123
 124   fPartEta->Scale(1.0/fEvents);
 125   fPartEta->Scale(1.0/fPartEta->GetBinWidth(1));
 126
 127   TCanvas* canvas = new TCanvas("control", "control", 900, 450);
 128   canvas->Divide(2, 1);
 129
 130   canvas->cd(1);
 131   fVertex->Draw();
 132
 133   canvas->cd(2);
 134   fPartEta->Draw();
 135 }