o) new way of calculating the final dndeta
[u/mrichter/AliRoot.git] / 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(), 1);
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 }