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