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