Added selector for dNdEtaAnalysis, use with testAnalysis2.C
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaAnalysisSelector.cxx
1 #include "AlidNdEtaAnalysisSelector.h"
2
3 #include <TStyle.h>
4 #include <TSystem.h>
5 #include <TCanvas.h>
6 #include <TParticle.h>
7 #include <TParticlePDG.h>
8 #include <TVector3.h>
9
10 #include <AliLog.h>
11 #include <AliGenEventHeader.h>
12 #include <AliTracker.h>
13
14 #include "../esdTrackCuts/AliESDtrackCuts.h"
15 #include "dNdEtaCorrection.h"
16 #include "dNdEtaAnalysis.h"
17
18 ClassImp(AlidNdEtaAnalysisSelector)
19
20 AlidNdEtaAnalysisSelector::AlidNdEtaAnalysisSelector(TTree *) :
21   AliSelector(),
22   fEsdTrackCuts(0),
23   fdNdEtaAnalysis(0),
24   fdNdEtaCorrection(0),
25   fdNdEtaAnalysisFinal(0)
26 {
27   //
28   // Constructor. Initialization of pointers
29   //
30 }
31
32 AlidNdEtaAnalysisSelector::~AlidNdEtaAnalysisSelector()
33 {
34   //
35   // Destructor
36   //
37
38   // histograms are in the output list and deleted when the output
39   // list is deleted by the TSelector dtor
40 }
41
42 void AlidNdEtaAnalysisSelector::SlaveBegin(TTree * tree)
43 {
44   // The SlaveBegin() function is called after the Begin() function.
45   // When running with PROOF SlaveBegin() is called on each slave server.
46   // The tree argument is deprecated (on PROOF 0 is passed).
47
48   AliSelector::SlaveBegin(tree);
49
50   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta");
51
52   if (fChain)
53   {
54     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fChain->GetUserInfo()->FindObject("AliESDtrackCuts"));
55     fdNdEtaCorrection = dynamic_cast<dNdEtaCorrection*> (fChain->GetUserInfo()->FindObject("dNdEtaCorrection"));
56   }
57
58   if (!fEsdTrackCuts)
59     printf("ERROR: Could not read EsdTrackCuts from user info\n");
60
61   if (!fEsdTrackCuts)
62     printf("ERROR: Could not read dNdEtaCorrection from user info\n");
63
64   AliLog::SetClassDebugLevel("AliESDtrackCuts", 1);
65 }
66
67 Bool_t AlidNdEtaAnalysisSelector::Process(Long64_t entry)
68 {
69   // The Process() function is called for each entry in the tree (or possibly
70   // keyed object in the case of PROOF) to be processed. The entry argument
71   // specifies which entry in the currently loaded tree is to be processed.
72   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
73   // to read either all or the required parts of the data. When processing
74   // keyed objects with PROOF, the object is already loaded and is available
75   // via the fObject pointer.
76   //
77   // This function should contain the "body" of the analysis. It can contain
78   // simple or elaborate selection criteria, run algorithms on the data
79   // of the event and typically fill histograms.
80
81   // WARNING when a selector is used with a TChain, you must use
82   //  the pointer to the current TTree to call GetEntry(entry).
83   //  The entry is always the local entry number in the current tree.
84   //  Assuming that fChain is the pointer to the TChain being processed,
85   //  use fChain->GetTree()->GetEntry(entry).
86
87   if (AliSelector::Process(entry) == kFALSE)
88     return kFALSE;
89
90   // Check prerequisites
91   if (!fESD || !fEsdTrackCuts)
92     return kFALSE;
93
94   // ########################################################
95   // get the EDS vertex
96   const AliESDVertex* vtxESD = fESD->GetVertex();
97
98   // the vertex should be reconstructed
99   if (strcmp(vtxESD->GetName(),"default")==0)
100     return kTRUE;
101
102   Double_t vtx_res[3];
103   vtx_res[0] = vtxESD->GetXRes();
104   vtx_res[1] = vtxESD->GetYRes();
105   vtx_res[2] = vtxESD->GetZRes();
106
107   // the resolution should be reasonable???
108   if (vtx_res[2]==0 || vtx_res[2]>0.1)
109     return kTRUE;
110
111   Double_t vtx[3];
112   vtxESD->GetXYZ(vtx);
113
114   // ########################################################
115   // loop over esd tracks
116   Int_t nTracks = fESD->GetNumberOfTracks();
117   for (Int_t t=0; t<nTracks; t++)
118   {
119     AliESDtrack* esdTrack = fESD->GetTrack(t);
120     if (!esdTrack)
121     {
122       printf("ERROR: Could not retrieve track %d.\n", t);
123       continue;
124     }
125
126     // cut the esd track?
127     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
128       continue;
129
130     Double_t p[3];
131     esdTrack->GetConstrainedPxPyPz(p); // ### TODO or GetInnerPxPyPy / GetOuter
132     TVector3 vector(p);
133
134     Float_t theta = vector.Theta();
135     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
136
137     Float_t correction = fdNdEtaCorrection->GetCorrection(vtx[2], eta);
138
139     fdNdEtaAnalysis->FillTrack(vtx[2], eta, correction);
140
141   } // end of track loop
142
143   // for event count per vertex
144   fdNdEtaAnalysis->FillEvent(vtx[2]);
145
146   return kTRUE;
147 }
148
149 void AlidNdEtaAnalysisSelector::SlaveTerminate()
150 {
151   // The SlaveTerminate() function is called after all entries or objects
152   // have been processed. When running with PROOF SlaveTerminate() is called
153   // on each slave server.
154
155   AliSelector::SlaveTerminate();
156
157   // Add the histograms to the output on each slave server
158   if (!fOutput)
159   {
160     printf("ERROR: Output list not initialized\n");
161     return;
162   }
163
164   fOutput->Add(fdNdEtaAnalysis->GetEtaVsVtxHistogram());
165   fOutput->Add(fdNdEtaAnalysis->GetEtaVsVtxUncorrectedHistogram());
166   fOutput->Add(fdNdEtaAnalysis->GetVtxHistogram());
167
168   fdNdEtaAnalysis->GetVtxHistogram()->Print();
169   fOutput->Print();
170 }
171
172 void AlidNdEtaAnalysisSelector::Terminate()
173 {
174   // The Terminate() function is the last function to be called during
175   // a query. It always runs on the client, it can be used to present
176   // the results graphically or save the results to file.
177
178   AliSelector::Terminate();
179
180   TH2F* etaVsVtxHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("eta_vs_vtx"));
181   TH2F* etaVsVtxUncorrectedHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("eta_vs_vtx_uncorrected"));
182   TH1D* vtxHistogram = dynamic_cast<TH1D*> (fOutput->FindObject("vtx"));
183
184   if (!etaVsVtxHistogram || !vtxHistogram || !etaVsVtxUncorrectedHistogram)
185   {
186     printf("ERROR: Histograms not available %p %p %p\n", (void*) etaVsVtxHistogram, (void*) etaVsVtxUncorrectedHistogram, (void*) vtxHistogram);
187     return;
188   }
189
190   fdNdEtaAnalysisFinal = new dNdEtaAnalysis("dNdEtaResult");
191
192   fdNdEtaAnalysisFinal->SetEtaVsVtxHistogram(etaVsVtxHistogram);
193   fdNdEtaAnalysisFinal->SetEtaVsVtxUncorrectedHistogram(etaVsVtxUncorrectedHistogram);
194   fdNdEtaAnalysisFinal->SetVtxHistogram(vtxHistogram);
195
196   fdNdEtaAnalysisFinal->Finish();
197
198   TFile* fout = new TFile("out.root","RECREATE");
199
200   fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
201   fdNdEtaCorrection->SaveHistograms();
202   fdNdEtaAnalysisFinal->SaveHistograms();
203
204   fout->Write();
205   fout->Close();
206
207   fdNdEtaAnalysisFinal->DrawHistograms();
208 }