changing from printf to AliDebug
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrectionSelector.cxx
CommitLineData
79ab56b9 1#include "AlidNdEtaCorrectionSelector.h"
539b6cb4 2
3#include <TStyle.h>
4#include <TSystem.h>
5#include <TCanvas.h>
6#include <TParticle.h>
7#include <TParticlePDG.h>
8
9#include <AliLog.h>
79ab56b9 10#include <AliGenEventHeader.h>
539b6cb4 11#include <AliTracker.h>
12
79ab56b9 13#include "../esdTrackCuts/AliESDtrackCuts.h"
14#include "dNdEtaCorrection.h"
539b6cb4 15
79ab56b9 16ClassImp(AlidNdEtaCorrectionSelector)
539b6cb4 17
79ab56b9 18AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector(TTree *) :
539b6cb4 19 AliSelector(),
20 fEsdTrackCuts(0),
75e130df 21 fdNdEtaCorrection(0),
22 fdNdEtaCorrectionFinal(0)
539b6cb4 23{
79ab56b9 24 //
539b6cb4 25 // Constructor. Initialization of pointers
79ab56b9 26 //
539b6cb4 27}
28
79ab56b9 29AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
539b6cb4 30{
79ab56b9 31 //
32 // Destructor
33 //
539b6cb4 34
35 // histograms are in the output list and deleted when the output
36 // list is deleted by the TSelector dtor
37}
38
79ab56b9 39void AlidNdEtaCorrectionSelector::Begin(TTree * tree)
539b6cb4 40{
41 // The Begin() function is called at the start of the query.
42 // When running with PROOF Begin() is only called on the client.
43 // The tree argument is deprecated (on PROOF 0 is passed).
44
45 AliSelector::Begin(tree);
539b6cb4 46}
47
79ab56b9 48void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
539b6cb4 49{
50 // The SlaveBegin() function is called after the Begin() function.
51 // When running with PROOF SlaveBegin() is called on each slave server.
52 // The tree argument is deprecated (on PROOF 0 is passed).
53
54 AliSelector::SlaveBegin(tree);
55
79ab56b9 56 fdNdEtaCorrection = new dNdEtaCorrection();
539b6cb4 57
79ab56b9 58 if (fChain)
59 fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fChain->GetUserInfo()->FindObject("AliESDtrackCuts"));
539b6cb4 60
79ab56b9 61 if (!fEsdTrackCuts)
b8e8168f 62 AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info");
539b6cb4 63}
64
79ab56b9 65Bool_t AlidNdEtaCorrectionSelector::Notify()
539b6cb4 66{
67 // The Notify() function is called when a new file is opened. This
68 // can be either for a new TTree in a TChain or when when a new TTree
69 // is started when using PROOF. Typically here the branch pointers
70 // will be retrieved. It is normaly not necessary to make changes
71 // to the generated code, but the routine can be extended by the
72 // user if needed.
73
74 if (AliSelector::Notify() == kFALSE)
75 return kFALSE;
76
539b6cb4 77 return kTRUE;
78}
79
79ab56b9 80Bool_t AlidNdEtaCorrectionSelector::IsPrimary(const TParticle* aParticle, Int_t aTotalPrimaries)
539b6cb4 81{
79ab56b9 82 //
83 // Returns if the given particle is a primary particle
84 // This function or a equivalent should be available in some common place of AliRoot
85 //
86
539b6cb4 87 // if the particle has a daughter primary, we do not want to count it
88 if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
89 return kFALSE;
90
91 Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
92
93 // skip quarks and gluon
94 if (pdgCode > 10 && pdgCode != 21)
95 return kTRUE;
96
97 return kFALSE;
98}
99
79ab56b9 100Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
539b6cb4 101{
102 // The Process() function is called for each entry in the tree (or possibly
103 // keyed object in the case of PROOF) to be processed. The entry argument
104 // specifies which entry in the currently loaded tree is to be processed.
105 // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
106 // to read either all or the required parts of the data. When processing
107 // keyed objects with PROOF, the object is already loaded and is available
108 // via the fObject pointer.
109 //
110 // This function should contain the "body" of the analysis. It can contain
111 // simple or elaborate selection criteria, run algorithms on the data
112 // of the event and typically fill histograms.
113
114 // WARNING when a selector is used with a TChain, you must use
115 // the pointer to the current TTree to call GetEntry(entry).
116 // The entry is always the local entry number in the current tree.
117 // Assuming that fChain is the pointer to the TChain being processed,
118 // use fChain->GetTree()->GetEntry(entry).
119
120 if (AliSelector::Process(entry) == kFALSE)
121 return kFALSE;
122
b8e8168f 123 // check prerequesites
124 if (!fESD)
125 {
126 AliDebug(AliLog::kError, "ESD branch not available");
127 return kFALSE;
128 }
129
130 if (!fHeader)
131 {
132 AliDebug(AliLog::kError, "Header branch not available");
133 return kFALSE;
134 }
135
136 if (!fEsdTrackCuts)
137 {
138 AliDebug(AliLog::kError, "fESDTrackCuts not available");
539b6cb4 139 return kFALSE;
b8e8168f 140 }
539b6cb4 141
142 // ########################################################
143 // get the EDS vertex
144 const AliESDVertex* vtxESD = fESD->GetVertex();
145
79ab56b9 146 // the vertex should be reconstructed
147 if (strcmp(vtxESD->GetName(),"default")==0)
148 return kTRUE;
539b6cb4 149
79ab56b9 150 Double_t vtx_res[3];
539b6cb4 151 vtx_res[0] = vtxESD->GetXRes();
152 vtx_res[1] = vtxESD->GetYRes();
153 vtx_res[2] = vtxESD->GetZRes();
154
539b6cb4 155 // the resolution should be reasonable???
156 if (vtx_res[2]==0 || vtx_res[2]>0.1)
157 return kTRUE;
158
159 // ########################################################
160 // get the MC vertex
79ab56b9 161 AliGenEventHeader* genHeader = fHeader->GenEventHeader();
539b6cb4 162
163 TArrayF vtxMC(3);
164 genHeader->PrimaryVertex(vtxMC);
539b6cb4 165
166 // ########################################################
167 // loop over mc particles
168 TTree* particleTree = GetKinematics();
169 TParticle* particle = 0;
170 particleTree->SetBranchAddress("Particles", &particle);
171
172 Int_t nPrim = fHeader->GetNprimary();
173 Int_t nTotal = fHeader->GetNtrack();
174
175 for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
176 {
177 particleTree->GetEntry(i_mc);
178
179 if (!particle)
180 continue;
181
182 if (strcmp(particle->GetName(),"XXX") == 0)
183 {
b8e8168f 184 AliDebug(AliLog::kWarning, Form("WARNING: There is a particle named XXX (%d).", i_mc));
539b6cb4 185 continue;
186 }
187
188 TParticlePDG* pdgPart = particle->GetPDG();
189
190 if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
191 {
b8e8168f 192 AliDebug(AliLog::kError, Form("WARNING: There is a particle with an unknown particle class (%d pdg code %d).", i_mc, particle->GetPdgCode()));
539b6cb4 193 continue;
194 }
195
196 if (IsPrimary(particle, nPrim) == kFALSE)
197 continue;
198
199 if (pdgPart->Charge() == 0)
200 continue;
201
79ab56b9 202 fdNdEtaCorrection->FillGene(vtxMC[2], particle->Eta());
539b6cb4 203
204 }// end of mc particle
205
206 // ########################################################
207 // loop over esd tracks
208 Int_t nTracks = fESD->GetNumberOfTracks();
209 for (Int_t t=0; t<nTracks; t++)
210 {
211 AliESDtrack* esdTrack = fESD->GetTrack(t);
212
213 // cut the esd track?
214 if (!fEsdTrackCuts->AcceptTrack(esdTrack))
215 continue;
216
539b6cb4 217 // using the eta of the mc particle
218 Int_t label = TMath::Abs(esdTrack->GetLabel());
79ab56b9 219 if (label == 0)
539b6cb4 220 {
b8e8168f 221 AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
539b6cb4 222 continue;
223 }
224 particleTree->GetEntry(nTotal - nPrim + label);
225
79ab56b9 226 fdNdEtaCorrection->FillMeas(vtxMC[2], particle->Eta());
539b6cb4 227
228 } // end of track loop
229
230 return kTRUE;
231}
232
79ab56b9 233void AlidNdEtaCorrectionSelector::SlaveTerminate()
539b6cb4 234{
235 // The SlaveTerminate() function is called after all entries or objects
236 // have been processed. When running with PROOF SlaveTerminate() is called
237 // on each slave server.
238
239 AliSelector::SlaveTerminate();
240
241 // Add the histograms to the output on each slave server
242 if (!fOutput)
243 {
b8e8168f 244 AliDebug(AliLog::kError, "ERROR: Output list not initialized");
539b6cb4 245 return;
246 }
79ab56b9 247
248 fOutput->Add(fdNdEtaCorrection->GetGeneratedHistogram());
249 fOutput->Add(fdNdEtaCorrection->GetMeasuredHistogram());
539b6cb4 250}
251
79ab56b9 252void AlidNdEtaCorrectionSelector::Terminate()
539b6cb4 253{
254 // The Terminate() function is the last function to be called during
255 // a query. It always runs on the client, it can be used to present
256 // the results graphically or save the results to file.
257
258 AliSelector::Terminate();
259
79ab56b9 260 fdNdEtaCorrectionFinal = new dNdEtaCorrection();
261 TH2F* measuredHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_meas"));
262 TH2F* generatedHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_gene"));
263 if (!measuredHistogram || !generatedHistogram)
264 {
b8e8168f 265 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) generatedHistogram, (void*) measuredHistogram));
79ab56b9 266 return;
267 }
268 fdNdEtaCorrectionFinal->SetGeneratedHistogram(generatedHistogram);
269 fdNdEtaCorrectionFinal->SetMeasuredHistogram(measuredHistogram);
270
271 fdNdEtaCorrectionFinal->Finish();
539b6cb4 272
273 TFile* fout = new TFile("correction_map.root","RECREATE");
274
275 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
79ab56b9 276 fdNdEtaCorrectionFinal->SaveHistograms();
539b6cb4 277
278 fout->Write();
279 fout->Close();
280
79ab56b9 281 fdNdEtaCorrectionFinal->DrawHistograms();
539b6cb4 282}