Added selector for dNdEtaAnalysis, use with testAnalysis2.C
[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)
62 printf("ERROR: Could not read EsdTrackCuts from user info\n");
539b6cb4 63
64 AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
65}
66
79ab56b9 67Bool_t AlidNdEtaCorrectionSelector::Notify()
539b6cb4 68{
69 // The Notify() function is called when a new file is opened. This
70 // can be either for a new TTree in a TChain or when when a new TTree
71 // is started when using PROOF. Typically here the branch pointers
72 // will be retrieved. It is normaly not necessary to make changes
73 // to the generated code, but the routine can be extended by the
74 // user if needed.
75
76 if (AliSelector::Notify() == kFALSE)
77 return kFALSE;
78
539b6cb4 79 return kTRUE;
80}
81
79ab56b9 82Bool_t AlidNdEtaCorrectionSelector::IsPrimary(const TParticle* aParticle, Int_t aTotalPrimaries)
539b6cb4 83{
79ab56b9 84 //
85 // Returns if the given particle is a primary particle
86 // This function or a equivalent should be available in some common place of AliRoot
87 //
88
539b6cb4 89 // if the particle has a daughter primary, we do not want to count it
90 if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
91 return kFALSE;
92
93 Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
94
95 // skip quarks and gluon
96 if (pdgCode > 10 && pdgCode != 21)
97 return kTRUE;
98
99 return kFALSE;
100}
101
79ab56b9 102Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
539b6cb4 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 fChain is the pointer to the TChain being processed,
120 // use fChain->GetTree()->GetEntry(entry).
121
122 if (AliSelector::Process(entry) == kFALSE)
123 return kFALSE;
124
125 if (!fESD || !fHeader)
126 return kFALSE;
127
128 // ########################################################
129 // get the EDS vertex
130 const AliESDVertex* vtxESD = fESD->GetVertex();
131
79ab56b9 132 // the vertex should be reconstructed
133 if (strcmp(vtxESD->GetName(),"default")==0)
134 return kTRUE;
539b6cb4 135
79ab56b9 136 Double_t vtx_res[3];
539b6cb4 137 vtx_res[0] = vtxESD->GetXRes();
138 vtx_res[1] = vtxESD->GetYRes();
139 vtx_res[2] = vtxESD->GetZRes();
140
539b6cb4 141 // the resolution should be reasonable???
142 if (vtx_res[2]==0 || vtx_res[2]>0.1)
143 return kTRUE;
144
145 // ########################################################
146 // get the MC vertex
79ab56b9 147 AliGenEventHeader* genHeader = fHeader->GenEventHeader();
539b6cb4 148
149 TArrayF vtxMC(3);
150 genHeader->PrimaryVertex(vtxMC);
539b6cb4 151
152 // ########################################################
153 // loop over mc particles
154 TTree* particleTree = GetKinematics();
155 TParticle* particle = 0;
156 particleTree->SetBranchAddress("Particles", &particle);
157
158 Int_t nPrim = fHeader->GetNprimary();
159 Int_t nTotal = fHeader->GetNtrack();
160
161 for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
162 {
163 particleTree->GetEntry(i_mc);
164
165 if (!particle)
166 continue;
167
168 if (strcmp(particle->GetName(),"XXX") == 0)
169 {
170 printf("WARNING: There is a particle named XXX (%d).\n", i_mc);
171 continue;
172 }
173
174 TParticlePDG* pdgPart = particle->GetPDG();
175
176 if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
177 {
178 printf("WARNING: There is a particle with an unknown particle class (%d pdg code %d).\n", i_mc, particle->GetPdgCode());
179 continue;
180 }
181
182 if (IsPrimary(particle, nPrim) == kFALSE)
183 continue;
184
185 if (pdgPart->Charge() == 0)
186 continue;
187
79ab56b9 188 fdNdEtaCorrection->FillGene(vtxMC[2], particle->Eta());
539b6cb4 189
190 }// end of mc particle
191
192 // ########################################################
193 // loop over esd tracks
194 Int_t nTracks = fESD->GetNumberOfTracks();
195 for (Int_t t=0; t<nTracks; t++)
196 {
197 AliESDtrack* esdTrack = fESD->GetTrack(t);
198
199 // cut the esd track?
200 if (!fEsdTrackCuts->AcceptTrack(esdTrack))
201 continue;
202
539b6cb4 203 // using the eta of the mc particle
204 Int_t label = TMath::Abs(esdTrack->GetLabel());
79ab56b9 205 if (label == 0)
539b6cb4 206 {
79ab56b9 207 printf("WARNING: cannot find corresponding mc part for track %d.", t);
539b6cb4 208 continue;
209 }
210 particleTree->GetEntry(nTotal - nPrim + label);
211
79ab56b9 212 fdNdEtaCorrection->FillMeas(vtxMC[2], particle->Eta());
539b6cb4 213
214 } // end of track loop
215
216 return kTRUE;
217}
218
79ab56b9 219void AlidNdEtaCorrectionSelector::SlaveTerminate()
539b6cb4 220{
221 // The SlaveTerminate() function is called after all entries or objects
222 // have been processed. When running with PROOF SlaveTerminate() is called
223 // on each slave server.
224
225 AliSelector::SlaveTerminate();
226
227 // Add the histograms to the output on each slave server
228 if (!fOutput)
229 {
230 printf("ERROR: Output list not initialized\n");
231 return;
232 }
79ab56b9 233
234 fOutput->Add(fdNdEtaCorrection->GetGeneratedHistogram());
235 fOutput->Add(fdNdEtaCorrection->GetMeasuredHistogram());
539b6cb4 236}
237
79ab56b9 238void AlidNdEtaCorrectionSelector::Terminate()
539b6cb4 239{
240 // The Terminate() function is the last function to be called during
241 // a query. It always runs on the client, it can be used to present
242 // the results graphically or save the results to file.
243
244 AliSelector::Terminate();
245
79ab56b9 246 fdNdEtaCorrectionFinal = new dNdEtaCorrection();
247 TH2F* measuredHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_meas"));
248 TH2F* generatedHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_gene"));
249 if (!measuredHistogram || !generatedHistogram)
250 {
251 printf("ERROR: Histograms not available %p %p\n", (void*) generatedHistogram, (void*) measuredHistogram);
252 return;
253 }
254 fdNdEtaCorrectionFinal->SetGeneratedHistogram(generatedHistogram);
255 fdNdEtaCorrectionFinal->SetMeasuredHistogram(measuredHistogram);
256
257 fdNdEtaCorrectionFinal->Finish();
539b6cb4 258
259 TFile* fout = new TFile("correction_map.root","RECREATE");
260
261 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
79ab56b9 262 fdNdEtaCorrectionFinal->SaveHistograms();
539b6cb4 263
264 fout->Write();
265 fout->Close();
266
79ab56b9 267 fdNdEtaCorrectionFinal->DrawHistograms();
539b6cb4 268}