1) Added "Makefile support": LinkDef + pkg file, made header files compliant
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaEffSelector.cxx
1 #include "AlidNdEtaEffSelector.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 <../PYTHIA6/AliGenPythiaEventHeader.h>
11 #include <../TPC/AliTPCtrack.h>
12 #include <AliTracker.h>
13
14 #include <iostream>
15 using namespace std;
16
17 ClassImp(AlidNdEtaEffSelector)
18
19 AlidNdEtaEffSelector::AlidNdEtaEffSelector(TTree *) :
20   AliSelector(),
21   fEsdTrackCuts(0),
22   fdNdEtaCorrection(0)
23 {
24   // Constructor. Initialization of pointers
25 }
26
27 AlidNdEtaEffSelector::~AlidNdEtaEffSelector()
28 {
29   // Remove all pointers
30
31   // histograms are in the output list and deleted when the output
32   // list is deleted by the TSelector dtor
33 }
34
35 void AlidNdEtaEffSelector::Begin(TTree * tree)
36 {
37   // The Begin() function is called at the start of the query.
38   // When running with PROOF Begin() is only called on the client.
39   // The tree argument is deprecated (on PROOF 0 is passed).
40
41   AliSelector::Begin(tree);
42
43   fdNdEtaCorrection = new dNdEtaCorrection();
44 }
45
46 void AlidNdEtaEffSelector::SlaveBegin(TTree * tree)
47 {
48   // The SlaveBegin() function is called after the Begin() function.
49   // When running with PROOF SlaveBegin() is called on each slave server.
50   // The tree argument is deprecated (on PROOF 0 is passed).
51
52   AliSelector::SlaveBegin(tree);
53
54   fEsdTrackCuts = new ESDtrackQualityCuts();
55   fEsdTrackCuts->DefineHistograms(1);
56
57   fEsdTrackCuts->SetMinNClustersTPC(50);
58   fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
59   fEsdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
60   fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
61
62   fEsdTrackCuts->SetMinNsigmaToVertex(3);
63   fEsdTrackCuts->SetAcceptKingDaughters(kFALSE);
64
65   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
66 }
67
68 Bool_t AlidNdEtaEffSelector::Notify()
69 {
70   // The Notify() function is called when a new file is opened. This
71   // can be either for a new TTree in a TChain or when when a new TTree
72   // is started when using PROOF. Typically here the branch pointers
73   // will be retrieved. It is normaly not necessary to make changes
74   // to the generated code, but the routine can be extended by the
75   // user if needed.
76
77   if (AliSelector::Notify() == kFALSE)
78     return kFALSE;
79
80   // ########################################################
81   // Magnetic field
82   AliTracker::SetFieldMap(GetAliRun()->Field(), kTRUE); // kTRUE means uniform magnetic field
83
84   return kTRUE;
85 }
86
87 Bool_t AlidNdEtaEffSelector::IsPrimary(const TParticle* aParticle, Int_t aTotalPrimaries)
88 {
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
102 Bool_t AlidNdEtaEffSelector::Process(Long64_t entry)
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
132   Double_t vtx[3];
133   Double_t vtx_res[3];
134   vtxESD->GetXYZ(vtx);
135
136   vtx_res[0] = vtxESD->GetXRes();
137   vtx_res[1] = vtxESD->GetYRes();
138   vtx_res[2] = vtxESD->GetZRes();
139
140   // the vertex should be reconstructed
141   if (strcmp(vtxESD->GetName(),"default")==0)
142     return kTRUE;
143
144   // the resolution should be reasonable???
145   if (vtx_res[2]==0 || vtx_res[2]>0.1)
146     return kTRUE;
147
148   // ########################################################
149   // get the MC vertex
150   AliGenPythiaEventHeader* genHeader = (AliGenPythiaEventHeader*) fHeader->GenEventHeader();
151
152   TArrayF vtxMC(3);
153   genHeader->PrimaryVertex(vtxMC);
154   vtx[0] = vtxMC[0];
155   vtx[1] = vtxMC[1];
156   vtx[2] = vtxMC[2];
157
158   // ########################################################
159   // loop over mc particles
160   TTree* particleTree = GetKinematics();
161   TParticle* particle = 0;
162   particleTree->SetBranchAddress("Particles", &particle);
163
164   Int_t nPrim  = fHeader->GetNprimary();
165   Int_t nTotal = fHeader->GetNtrack();
166
167   for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
168   {
169     particleTree->GetEntry(i_mc);
170
171     if (!particle)
172       continue;
173
174     if (strcmp(particle->GetName(),"XXX") == 0)
175     {
176       printf("WARNING: There is a particle named XXX (%d).\n", i_mc);
177       continue;
178     }
179
180     TParticlePDG* pdgPart = particle->GetPDG();
181
182     if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
183     {
184       printf("WARNING: There is a particle with an unknown particle class (%d pdg code %d).\n", i_mc, particle->GetPdgCode());
185       continue;
186     }
187
188     if (IsPrimary(particle, nPrim) == kFALSE)
189       continue;
190
191     if (pdgPart->Charge() == 0)
192       continue;
193
194     fdNdEtaCorrection->FillGene(vtx[2], particle->Eta());
195
196   }// end of mc particle
197
198   // ########################################################
199   // loop over esd tracks
200   Int_t nTracks = fESD->GetNumberOfTracks();
201   for (Int_t t=0; t<nTracks; t++)
202   {
203     AliESDtrack* esdTrack = fESD->GetTrack(t);
204
205     // cut the esd track?
206     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
207       continue;
208
209     AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
210     if (tpcTrack->GetAlpha()==0)
211     {
212       cout << " Warning esd track alpha = 0" << endl;
213       continue;
214     }
215
216     //Float_t theta = tpcTrack->Theta();
217     //Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
218
219     // using the eta of the mc particle
220     Int_t label = TMath::Abs(esdTrack->GetLabel());
221     if (label<0) 
222     {
223       cout << " cannot find corresponding mc part !!! " << label << endl;
224       continue;
225     }
226     particleTree->GetEntry(nTotal - nPrim + label);
227
228     fdNdEtaCorrection->FillMeas(vtx[2], particle->Eta());
229
230   } // end of track loop
231
232   return kTRUE;
233 }
234
235 void AlidNdEtaEffSelector::SlaveTerminate()
236 {
237   // The SlaveTerminate() function is called after all entries or objects
238   // have been processed. When running with PROOF SlaveTerminate() is called
239   // on each slave server.
240
241   AliSelector::SlaveTerminate();
242
243   // Add the histograms to the output on each slave server
244   if (!fOutput)
245   {
246     printf("ERROR: Output list not initialized\n");
247     return;
248   }
249 }
250
251 void AlidNdEtaEffSelector::Terminate()
252 {
253   // The Terminate() function is the last function to be called during
254   // a query. It always runs on the client, it can be used to present
255   // the results graphically or save the results to file.
256
257   AliSelector::Terminate();
258
259   fdNdEtaCorrection->Finish();
260
261   TFile* fout = new TFile("correction_map.root","RECREATE");
262
263   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
264   fdNdEtaCorrection->SaveHistograms();
265
266   fout->Write();
267   fout->Close();
268
269   fdNdEtaCorrection->DrawHistograms();
270 }