adding multipicity from mc, and creation of correlation map for the corrections
[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 <TH1F.h>
12 #include <TH2F.h>
13 #include <TParticle.h>
14
15 #include <AliLog.h>
16 #include <AliESD.h>
17 #include <AliStack.h>
18
19 #include "esdTrackCuts/AliESDtrackCuts.h"
20 #include "AliPWG0Helper.h"
21
22 ClassImp(AliMultiplicityMCSelector)
23
24 AliMultiplicityMCSelector::AliMultiplicityMCSelector() :
25   AliSelectorRL(),
26   fMultiplicityESD(0),
27   fMultiplicityMC(0),
28   fCorrelation(0),
29   fEsdTrackCuts(0)
30 {
31   //
32   // Constructor. Initialization of pointers
33   //
34 }
35
36 AliMultiplicityMCSelector::~AliMultiplicityMCSelector()
37 {
38   //
39   // Destructor
40   //
41
42   // histograms are in the output list and deleted when the output
43   // list is deleted by the TSelector dtor
44 }
45
46 void AliMultiplicityMCSelector::Begin(TTree* tree)
47 {
48   // Begin function
49
50   AliSelectorRL::Begin(tree);
51
52   ReadUserObjects(tree);
53 }
54
55 void AliMultiplicityMCSelector::ReadUserObjects(TTree* tree)
56 {
57   // read the user objects, called from slavebegin and begin
58
59   if (!fEsdTrackCuts && fInput)
60     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
61
62   if (!fEsdTrackCuts && tree)
63     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
64
65   if (!fEsdTrackCuts)
66      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
67 }
68
69 void AliMultiplicityMCSelector::SlaveBegin(TTree* tree)
70 {
71   // The SlaveBegin() function is called after the Begin() function.
72   // When running with PROOF SlaveBegin() is called on each slave server.
73   // The tree argument is deprecated (on PROOF 0 is passed).
74
75   AliSelector::SlaveBegin(tree);
76
77   ReadUserObjects(tree);
78
79   fMultiplicityESD = new TH1F("fMultiplicityESD", "fMultiplicityESD;Ntracks;Count", 201, -0.5, 200.5);
80   fMultiplicityMC = new TH1F("fMultiplicityMC", "fMultiplicityMC;Npart;Count", 201, -0.5, 200.5);
81
82   fCorrelation = new TH2F("fCorrelation", "fCorrelation;Npart;Ntracks", 201, -0.5, 200.5, 201, -0.5, 200.5);
83 }
84
85 Bool_t AliMultiplicityMCSelector::Process(Long64_t entry)
86 {
87   // The Process() function is called for each entry in the tree (or possibly
88   // keyed object in the case of PROOF) to be processed. The entry argument
89   // specifies which entry in the currently loaded tree is to be processed.
90   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
91   // to read either all or the required parts of the data. When processing
92   // keyed objects with PROOF, the object is already loaded and is available
93   // via the fObject pointer.
94   //
95   // This function should contain the "body" of the analysis. It can contain
96   // simple or elaborate selection criteria, run algorithms on the data
97   // of the event and typically fill histograms.
98
99   // WARNING when a selector is used with a TChain, you must use
100   //  the pointer to the current TTree to call GetEntry(entry).
101   //  The entry is always the local entry number in the current tree.
102   //  Assuming that fTree is the pointer to the TChain being processed,
103   //  use fTree->GetTree()->GetEntry(entry).
104
105   if (AliSelectorRL::Process(entry) == kFALSE)
106     return kFALSE;
107
108   // Check prerequisites
109   if (!fESD)
110   {
111     AliDebug(AliLog::kError, "ESD branch not available");
112     return kFALSE;
113   }
114
115   if (!fEsdTrackCuts)
116   {
117     AliDebug(AliLog::kError, "fESDTrackCuts not available");
118     return kFALSE;
119   }
120
121   AliStack* stack = GetStack();
122   if (!stack)
123   {
124     AliDebug(AliLog::kError, "Stack not available");
125     return kFALSE;
126   }
127
128   if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
129     return kTRUE;
130
131   if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
132     return kTRUE;
133
134   // TODO put pt cut
135
136   // get number of "good" tracks from ESD
137   Int_t nESDTracks = fEsdTrackCuts->CountAcceptedTracks(fESD);
138   fMultiplicityESD->Fill(nESDTracks);
139
140   // get number of tracks from MC
141   // loop over mc particles
142   Int_t nPrim  = stack->GetNprimary();
143   Int_t nMCTracks = 0;
144   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
145   {
146     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
147
148     TParticle* particle = stack->Particle(iMc);
149
150     if (!particle)
151     {
152       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
153       continue;
154     }
155
156     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
157       continue;
158
159     Float_t eta = particle->Eta();
160
161     if (TMath::Abs(eta) > 0.9)
162       continue;
163
164     nMCTracks++;
165   }// end of mc particle
166
167   fMultiplicityMC->Fill(nMCTracks);
168
169   fCorrelation->Fill(nMCTracks, nESDTracks);
170
171   return kTRUE;
172 }
173
174 void AliMultiplicityMCSelector::SlaveTerminate()
175 {
176   // The SlaveTerminate() function is called after all entries or objects
177   // have been processed. When running with PROOF SlaveTerminate() is called
178   // on each slave server.
179
180   AliSelector::SlaveTerminate();
181
182   // Add the histograms to the output on each slave server
183   if (!fOutput)
184   {
185     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
186     return;
187   }
188
189   fOutput->Add(fMultiplicityESD);
190   fOutput->Add(fMultiplicityMC);
191   fOutput->Add(fCorrelation);
192 }
193
194 void AliMultiplicityMCSelector::Terminate()
195 {
196   // The Terminate() function is the last function to be called during
197   // a query. It always runs on the client, it can be used to present
198   // the results graphically or save the results to file.
199
200   AliSelector::Terminate();
201
202   fOutput->Print();
203
204   fMultiplicityESD = dynamic_cast<TH1F*> (fOutput->FindObject("fMultiplicityESD"));
205   fMultiplicityMC = dynamic_cast<TH1F*> (fOutput->FindObject("fMultiplicityMC"));
206   fCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fCorrelation"));
207
208   if (!fMultiplicityESD || !fMultiplicityMC || !fCorrelation)
209   {
210     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p %p", (void*) fMultiplicityESD, (void*) fMultiplicityMC, (void*) fCorrelation));
211     return;
212   }
213
214   TCanvas* canvas = new TCanvas("AliMultiplicityMCSelector", "AliMultiplicityMCSelector", 1000, 1000);
215   canvas->Divide(2, 2);
216
217   canvas->cd(1);
218   fMultiplicityESD->Draw();
219   fMultiplicityMC->SetLineColor(2);
220   fMultiplicityMC->Draw("SAME");
221
222   canvas->cd(2);
223   TH1F* ratio = dynamic_cast<TH1F*> (fMultiplicityESD->Clone("multiplicity_ratio"));
224   ratio->SetTitle("ratio;Ntracks;Nreco/Ngene");
225   ratio->Divide(fMultiplicityMC);
226   ratio->Draw();
227
228   canvas->cd(3);
229   fCorrelation->Draw("COLZ");
230
231   TFile* file = TFile::Open("multiplicity.root", "RECREATE");
232
233   fMultiplicityMC->Write();
234   fMultiplicityESD->Write();
235   fCorrelation->Write();
236
237   file->Close();
238 }