introducing multiplicity measurement with the ITS or TPC
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AliMultiplicityMCSelector.cxx
CommitLineData
43a9a462 1/* $Id$ */
2
3#include "AliMultiplicityMCSelector.h"
4
5#include <TStyle.h>
6#include <TSystem.h>
7#include <TCanvas.h>
8#include <TVector3.h>
9#include <TChain.h>
10#include <TFile.h>
43a9a462 11#include <TH2F.h>
0a173978 12#include <TH3F.h>
43a9a462 13#include <TParticle.h>
14
15#include <AliLog.h>
16#include <AliESD.h>
17#include <AliStack.h>
0a173978 18#include <AliHeader.h>
19#include <AliGenEventHeader.h>
20#include <AliMultiplicity.h>
43a9a462 21
22#include "esdTrackCuts/AliESDtrackCuts.h"
23#include "AliPWG0Helper.h"
0a173978 24#include "dNdEta/AliMultiplicityCorrection.h"
43a9a462 25
26ClassImp(AliMultiplicityMCSelector)
27
28AliMultiplicityMCSelector::AliMultiplicityMCSelector() :
29 AliSelectorRL(),
0a173978 30 fMultiplicity(0),
43a9a462 31 fEsdTrackCuts(0)
32{
33 //
34 // Constructor. Initialization of pointers
35 //
36}
37
38AliMultiplicityMCSelector::~AliMultiplicityMCSelector()
39{
40 //
41 // Destructor
42 //
43
44 // histograms are in the output list and deleted when the output
45 // list is deleted by the TSelector dtor
46}
47
48void AliMultiplicityMCSelector::Begin(TTree* tree)
49{
50 // Begin function
51
52 AliSelectorRL::Begin(tree);
53
54 ReadUserObjects(tree);
55}
56
57void AliMultiplicityMCSelector::ReadUserObjects(TTree* tree)
58{
59 // read the user objects, called from slavebegin and begin
60
61 if (!fEsdTrackCuts && fInput)
62 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
63
64 if (!fEsdTrackCuts && tree)
65 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
66
67 if (!fEsdTrackCuts)
68 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
69}
70
71void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
72{
73 // The SlaveBegin() function is called after the Begin() function.
74 // When running with PROOF SlaveBegin() is called on each slave server.
75 // The tree argument is deprecated (on PROOF 0 is passed).
76
77 AliSelector::SlaveBegin(tree);
78
79 ReadUserObjects(tree);
80
0a173978 81 fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
82}
83
84void AliMultiplicityMCSelector::Init(TTree* tree)
85{
86 // read the user objects
87
88 AliSelectorRL::Init(tree);
89
90 // TODO Enable only the needed branches
91 /*if (tree)
92 {
93 tree->SetBranchStatus("*", 0);
94 tree->SetBranchStatus("fTriggerMask", 1);
95 tree->SetBranchStatus("fSPDVertex*", 1);
96 tree->SetBranchStatus("fTracks.fLabel", 1);
43a9a462 97
0a173978 98 AliESDtrackCuts::EnableNeededBranches(tree);
99 }*/
43a9a462 100}
101
102Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
103{
104 // The Process() function is called for each entry in the tree (or possibly
105 // keyed object in the case of PROOF) to be processed. The entry argument
106 // specifies which entry in the currently loaded tree is to be processed.
107 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
108 // to read either all or the required parts of the data. When processing
109 // keyed objects with PROOF, the object is already loaded and is available
110 // via the fObject pointer.
111 //
112 // This function should contain the "body" of the analysis. It can contain
113 // simple or elaborate selection criteria, run algorithms on the data
114 // of the event and typically fill histograms.
115
116 // WARNING when a selector is used with a TChain, you must use
117 // the pointer to the current TTree to call GetEntry(entry).
118 // The entry is always the local entry number in the current tree.
119 // Assuming that fTree is the pointer to the TChain being processed,
120 // use fTree->GetTree()->GetEntry(entry).
121
122 if (AliSelectorRL::Process(entry) == kFALSE)
123 return kFALSE;
124
125 // Check prerequisites
126 if (!fESD)
127 {
128 AliDebug(AliLog::kError, "ESD branch not available");
129 return kFALSE;
130 }
131
132 if (!fEsdTrackCuts)
133 {
134 AliDebug(AliLog::kError, "fESDTrackCuts not available");
135 return kFALSE;
136 }
137
0a173978 138 AliHeader* header = GetHeader();
139 if (!header)
140 {
141 AliDebug(AliLog::kError, "Header not available");
142 return kFALSE;
143 }
144
43a9a462 145 AliStack* stack = GetStack();
146 if (!stack)
147 {
148 AliDebug(AliLog::kError, "Stack not available");
149 return kFALSE;
150 }
151
152 if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
153 return kTRUE;
154
155 if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
156 return kTRUE;
157
0a173978 158 // get the ESD vertex
159 const AliESDVertex* vtxESD = fESD->GetVertex();
160 Double_t vtx[3];
161 vtxESD->GetXYZ(vtx);
162
163 // get the MC vertex
164 AliGenEventHeader* genHeader = header->GenEventHeader();
165 TArrayF vtxMC(3);
166 genHeader->PrimaryVertex(vtxMC);
167
168 // get tracklets
169 const AliMultiplicity* mult = fESD->GetMultiplicity();
170 if (!mult)
171 {
172 AliDebug(AliLog::kError, "AliMultiplicity not available");
173 return kFALSE;
174 }
43a9a462 175
0a173978 176 //const Float_t kPtCut = 0.3;
43a9a462 177
178 // get number of tracks from MC
179 // loop over mc particles
180 Int_t nPrim = stack->GetNprimary();
0a173978 181
182 // tracks in different eta ranges
183 Int_t nMCTracks05 = 0;
184 Int_t nMCTracks10 = 0;
185 Int_t nMCTracks15 = 0;
186 Int_t nMCTracks20 = 0;
187 Int_t nMCTracksAll = 0;
188
43a9a462 189 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
190 {
191 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
192
193 TParticle* particle = stack->Particle(iMc);
194
195 if (!particle)
196 {
197 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
198 continue;
199 }
200
201 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
202 continue;
203
0a173978 204 //if (particle->Pt() < kPtCut)
205 // continue;
43a9a462 206
0a173978 207 if (TMath::Abs(particle->Eta()) < 0.5)
208 nMCTracks05++;
43a9a462 209
0a173978 210 if (TMath::Abs(particle->Eta()) < 1.0)
211 nMCTracks10++;
212
213 if (TMath::Abs(particle->Eta()) < 1.5)
214 nMCTracks15++;
215
216 if (TMath::Abs(particle->Eta()) < 2.0)
217 nMCTracks20++;
218
219 nMCTracksAll++;
43a9a462 220 }// end of mc particle
221
0a173978 222 // FAKE: put back vtxMC[2]
223 fMultiplicity->FillGenerated(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll);
224
225 Int_t nESDTracks05 = 0;
226 Int_t nESDTracks10 = 0;
227 Int_t nESDTracks15 = 0;
228 Int_t nESDTracks20 = 0;
229
230 // get multiplicity from ITS tracklets
231 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
232 {
233 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
234
235 // this removes non-tracklets. Very bad solution. SPD guys are working on better solution...
236 if (mult->GetDeltaPhi(i) < -1000)
237 continue;
238
239 Float_t theta = mult->GetTheta(i);
240 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
241
242 if (TMath::Abs(eta) < 0.5)
243 nESDTracks05++;
244
245 if (TMath::Abs(eta) < 1.0)
246 nESDTracks10++;
247
248 if (TMath::Abs(eta) < 1.5)
249 nESDTracks15++;
250
251 if (TMath::Abs(eta) < 2.0)
252 nESDTracks20++;
253 }
254
255 // get multiplicity from ESD tracks
256 /*TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
257 Int_t nGoodTracks = list->GetEntries();
258 // loop over esd tracks
259 for (Int_t i=0; i<nGoodTracks; i++)
260 {
261 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
262 if (!esdTrack)
263 {
264 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
265 continue;
266 }
267
268 Double_t p[3];
269 esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist
270 TVector3 vector(p);
271
272 Float_t theta = vector.Theta();
273 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
274 Float_t pt = vector.Pt();
275
276 if (pt < kPtCut)
277 continue;
278
279 if (TMath::Abs(eta) < 0.5)
280 nESDTracks05++;
43a9a462 281
0a173978 282 if (TMath::Abs(eta) < 1.0)
283 nESDTracks10++;
284
285 if (TMath::Abs(eta) < 1.5)
286 nESDTracks15++;
287
288 if (TMath::Abs(eta) < 2.0)
289 nESDTracks20++;
290 }*/
291
292 fMultiplicity->FillMeasured(vtxMC[2], nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
293
294 // TODO should this be vtx[2] or vtxMC[2] ?
295 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks15, nMCTracks20, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks15, nESDTracks20);
43a9a462 296
297 return kTRUE;
298}
299
300void AliMultiplicityMCSelector::SlaveTerminate()
301{
302 // The SlaveTerminate() function is called after all entries or objects
303 // have been processed. When running with PROOF SlaveTerminate() is called
304 // on each slave server.
305
306 AliSelector::SlaveTerminate();
307
308 // Add the histograms to the output on each slave server
309 if (!fOutput)
310 {
311 AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
312 return;
313 }
314
0a173978 315 fOutput->Add(fMultiplicity);
43a9a462 316}
317
318void AliMultiplicityMCSelector::Terminate()
319{
320 // The Terminate() function is the last function to be called during
321 // a query. It always runs on the client, it can be used to present
322 // the results graphically or save the results to file.
323
324 AliSelector::Terminate();
325
0a173978 326 fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
43a9a462 327
0a173978 328 if (!fMultiplicity)
43a9a462 329 {
0a173978 330 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
43a9a462 331 return;
332 }
333
f31d5d49 334 TFile* file = TFile::Open("multiplicityMC.root", "RECREATE");
43a9a462 335
0a173978 336 fMultiplicity->SaveHistograms();
43a9a462 337
338 file->Close();
0a173978 339
340 fMultiplicity->DrawHistograms();
43a9a462 341}