new Muon-Hadron correlation task in DHC (Constantin Loizides <cloizides@lbl.gov>...
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliMuonEffMC.cxx
CommitLineData
7e990b20 1// MUON track QA referring AliMuonEffMC.cxx
2// Author : Saehanseul Oh
3
4#include "AliMuonEffMC.h"
5
6#include <TList.h>
7#include <TH1D.h>
8#include <TH2F.h>
9#include <TH3F.h>
10#include <THn.h>
11#include <TChain.h>
12#include <TFile.h>
13
14#include "AliAnalysisManager.h"
15#include "AliAnalysisTask.h"
16#include "AliESDVertex.h"
17#include "AliESDEvent.h"
18#include "AliAODEvent.h"
19#include "AliESDMuonTrack.h"
20#include "AliESDVertex.h"
21#include "AliAODVertex.h"
22#include "AliCentrality.h"
23#include "AliVParticle.h"
24#include "AliMCParticle.h"
25#include "AliAnalysisHelperJetTasks.h"
26#include "AliMCEvent.h"
27
28using std::cout;
29using std::endl;
30
31ClassImp(AliMuonEffMC)
32
33//________________________________________________________________________
34AliMuonEffMC::AliMuonEffMC() :
35 AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
36 fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
37 fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
38 fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
39{
40 // Constructor
41 DefineInput(0, TChain::Class());
42 DefineOutput(1, TList::Class());
43
44 for(Int_t i=0; i<4; i++)
45 {
46 fHMuMotherGenPt[i] = 0x0;
47 fHMuMotherRecPt[i] = 0x0;
48 fHMuMotherGenPhi[i] = 0x0;
49 fHMuMotherRecPhi[i] = 0x0;
50 fHMuMotherGenEta[i] = 0x0;
51 fHMuMotherRecEta[i] = 0x0;
52 fHMuDCA[i] = 0x0;
53 }
54}
55
56//________________________________________________________________________
57AliMuonEffMC::AliMuonEffMC(const char *name) :
58 AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
59 fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
60 fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
61 fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
62{
63 // Constructor
64 DefineInput(0, TChain::Class());
65 DefineOutput(1, TList::Class());
66
67 for(Int_t i=0; i<4; i++)
68 {
69 fHMuMotherGenPt[i] = 0x0;
70 fHMuMotherRecPt[i] = 0x0;
71 fHMuMotherGenPhi[i] = 0x0;
72 fHMuMotherRecPhi[i] = 0x0;
73 fHMuMotherGenEta[i] = 0x0;
74 fHMuMotherRecEta[i] = 0x0;
75 fHMuDCA[i] = 0x0;
76 }
77}
78
79//________________________________________________________________________
80AliMuonEffMC::~AliMuonEffMC()
81{
82 //Destructor
83}
84
85//________________________________________________________________________
86void AliMuonEffMC::UserCreateOutputObjects()
87{
88 // Create histograms
89 // Called once (per slave on PROOF!)
90
91 fOutputList = new TList();
92 fOutputList->SetOwner(1);
93
94 fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
95 fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
96 fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
97 fHEventStat->GetXaxis()->SetBinLabel(3,"File");
98 fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt"); //!Global Trigger Single plus High p_T
99 fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt"); //!Global Trigger Single plus All p_T
100 fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt"); //!Global Trigger Single minus Low p_T
101 fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt"); //!Global Trigger Single minus High p_T
102 fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt"); //!Global Trigger Single minus All p_T
103 fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt"); //!Global Trigger Single undefined Low p_T
104 fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
105 fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt"); //!Global Trigger Single undefined All p_T
106 fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt"); //!Global Trigger UnlikeSign pair Low p_T
107 fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
108 fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt"); //!Global Trigger UnlikeSign pair All p_T
109 fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt"); //!Global Trigger LikeSign pair pair Low p_T
110 fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
111 fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt"); //!Global Trigger LikeSign pair pair All p_T
112 fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt"); //!Global Trigger Single plus Low p_T
113 fOutputList->Add(fHEventStat);
114
115 fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
116 fOutputList->Add(fHEvt);
117
118
119 Int_t iTrackBin[5];
120 Double_t* trackBins[5];
121 const char* trackAxisTitle[5];
122
123 // eta
124 Double_t etaBins[fNEtaBins+1];
125 for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
126 iTrackBin[0] = fNEtaBins;
127 trackBins[0] = etaBins;
128 trackAxisTitle[0] = "#eta";
129
130 // p_T
131 Double_t pTBins[fNpTBins+1];
132 for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
133 iTrackBin[1] = fNpTBins;
134 trackBins[1] = pTBins;
135 trackAxisTitle[1] = "p_{T} (GeV/c)";
136
137 // centrality
138 Double_t CentBins[fNCentBins+1];
139 for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
140 iTrackBin[2] = fNCentBins;
141 trackBins[2] = CentBins;
142 trackAxisTitle[2] = "Cent";
143
144 // Z-vertex
145 Double_t ZvtxBins[fNZvtxBins+1];
146 for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
147 iTrackBin[3] = fNZvtxBins;
148 trackBins[3] = ZvtxBins;
149 trackAxisTitle[3] = "Zvtx";
150
151 // phi
152 Double_t phiBins[fNPhiBins+1];
153 for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
154 iTrackBin[4] = fNPhiBins;
155 trackBins[4] = phiBins;
156 trackAxisTitle[4] = "#phi";
157
158 // THn for tracking efficiency
159 fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
160 for (Int_t i=0; i<5; i++)
161 {
162 fHMuonParGen->SetBinEdges(i, trackBins[i]);
163 fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
164 }
165 fHMuonParGen->Sumw2();
166 fOutputList->Add(fHMuonParGen);
167
168 fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
169 fHMuonDetGen->Sumw2();
170 fOutputList->Add(fHMuonDetGen);
171
172 fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
173 fHMuonDetRec->Sumw2();
174 fOutputList->Add(fHMuonDetRec);
175
176 fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
177 fHEtcDetRec->Sumw2();
178 fOutputList->Add(fHEtcDetRec);
179
180 if(fMDProcess)
181 {
182 const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
183
184 for(Int_t i=0; i<4; i++)
185 {
186 fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
187 fOutputList->Add(fHMuMotherGenPt[i]);
188
189 fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
190 fOutputList->Add(fHMuMotherRecPt[i]);
191
192 fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
193 fOutputList->Add(fHMuMotherGenPhi[i]);
194
195 fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
196 fOutputList->Add(fHMuMotherRecPhi[i]);
197
198 fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
199 fOutputList->Add(fHMuMotherGenEta[i]);
200
201 fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
202 fOutputList->Add(fHMuMotherRecEta[i]);
203
204 fHMuDCA[i] = new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
205 fOutputList->Add(fHMuDCA[i]);
206 }
207 }
208 PostData(1, fOutputList);
209}
210
211//________________________________________________________________________
212void AliMuonEffMC::UserExec(Option_t *)
213{
214 // Main loop, Called for each event
215 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
216 if (!fESD)
217 {
218 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
219 }
220
221 if (!fESD && !fAOD)
222 {
223 AliError("Neither fESD nor fAOD available");
224 return;
225 }
226
227 fHEventStat->Fill(0.5);
228 if(fIsMc)
229 {
230 fMC = MCEvent();
231 if (!fMC) {
232 printf("ERROR: fMC not available\n");
233 return;
234 }
235 }
236
237 const AliESDVertex* vertex = fESD->GetPrimaryVertex();
238 fZVertex = vertex->GetZ();
239 if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile("fCentralityEstimator");
240
241 // Fill Event histogram
242 fHEvt->Fill(fZVertex, fCentrality);
243
244 // Centrality, vertex, other event variables...
245 if (!VertexOk(fESD)) {
246 AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
247 return;
248 }
249
250 if (fCentrality > 100. || fCentrality < -1.5) {
251 AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
252 return;
253 }
254
255 if(fCentrality < 0) fCentrality = 1.0;
256 fHEventStat->Fill(1.5);
257
258 ULong64_t trigword=fESD->GetTriggerMask();
259
260 if (trigword & 0x01) fHEventStat->Fill(17.5);
261 if (trigword & 0x02) fHEventStat->Fill(3.5);
262 if (trigword & 0x04) fHEventStat->Fill(4.5);
263 if (trigword & 0x08) fHEventStat->Fill(5.5);
264 if (trigword & 0x010) fHEventStat->Fill(6.5);
265 if (trigword & 0x020) fHEventStat->Fill(7.5);
266 if (trigword & 0x040) fHEventStat->Fill(8.5);
267 if (trigword & 0x080) fHEventStat->Fill(9.5);
268 if (trigword & 0x100) fHEventStat->Fill(10.5);
269 if (trigword & 0x200) fHEventStat->Fill(11.5);
270 if (trigword & 0x400) fHEventStat->Fill(12.5);
271 if (trigword & 0x800) fHEventStat->Fill(13.5);
272 if (trigword & 0x1000) fHEventStat->Fill(14.5);
273 if (trigword & 0x2000) fHEventStat->Fill(15.5);
274 if (trigword & 0x4000) fHEventStat->Fill(16.5);
275
276 if(fIsMc)
277 {
278 for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
279 {
280 // loop on muon tracks
281 AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
282 if(muonTrack)
283 {
284 if(!IsGoodMUONtrack(*muonTrack)) continue;
285 Double_t fillArrayDetRec[5] = { muonTrack->Eta(), muonTrack->Pt(), fCentrality, fZVertex, muonTrack->Phi() };
286 fHMuonDetRec->Fill(fillArrayDetRec);
287
288 Int_t label = TMath::Abs(muonTrack->GetLabel());
289 AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(label);
290 if(TMath::Abs(McParticle->PdgCode()) != 13)
291 {
292 fHEtcDetRec->Fill(fillArrayDetRec);
293 continue;
294 }
295 Double_t fillArrayDetGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
296 fHMuonDetGen->Fill(fillArrayDetGen);
297
298 if(fMDProcess)
299 {
300 if(McParticle->GetMother() > 0)
301 {
302 AliMCParticle *MotherParticle = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
303 Int_t motherlabel = TMath::Abs(MotherParticle->PdgCode());
304
305 if(motherlabel==411 || motherlabel==413 || motherlabel==421 || motherlabel==423 || motherlabel==431 || motherlabel==433 || motherlabel==10413 || motherlabel==10411 || motherlabel==10423 || motherlabel==10421 || motherlabel==10433 || motherlabel==10431 || motherlabel==20413 || motherlabel==415 || motherlabel==20423 || motherlabel==425 || motherlabel==20433 || motherlabel==435)
306 {
307 fHMuMotherGenPt[2]->Fill(McParticle->Pt(), MotherParticle->Pt());
308 fHMuMotherRecPt[2]->Fill(muonTrack->Pt(), MotherParticle->Pt());
309 fHMuMotherGenPhi[2]->Fill(McParticle->Phi(), MotherParticle->Phi());
310 fHMuMotherRecPhi[2]->Fill(muonTrack->Phi(), MotherParticle->Phi());
311 fHMuMotherGenEta[2]->Fill(McParticle->Eta(), MotherParticle->Eta());
312 fHMuMotherRecEta[2]->Fill(muonTrack->Eta(), MotherParticle->Eta());
313 fHMuDCA[2]->Fill(muonTrack->GetDCA());
314 }
315 else if(motherlabel==211)
316 {
317 fHMuMotherGenPt[0]->Fill(McParticle->Pt(), MotherParticle->Pt());
318 fHMuMotherRecPt[0]->Fill(muonTrack->Pt(), MotherParticle->Pt());
319 fHMuMotherGenPhi[0]->Fill(McParticle->Phi(), MotherParticle->Phi());
320 fHMuMotherRecPhi[0]->Fill(muonTrack->Phi(), MotherParticle->Phi());
321 fHMuMotherGenEta[0]->Fill(McParticle->Eta(), MotherParticle->Eta());
322 fHMuMotherRecEta[0]->Fill(muonTrack->Eta(), MotherParticle->Eta());
323 fHMuDCA[0]->Fill(muonTrack->GetDCA());
324 }
325 else if(motherlabel==321)
326 {
327 fHMuMotherGenPt[1]->Fill(McParticle->Pt(), MotherParticle->Pt());
328 fHMuMotherRecPt[1]->Fill(muonTrack->Pt(), MotherParticle->Pt());
329 fHMuMotherGenPhi[1]->Fill(McParticle->Phi(), MotherParticle->Phi());
330 fHMuMotherRecPhi[1]->Fill(muonTrack->Phi(), MotherParticle->Phi());
331 fHMuMotherGenEta[1]->Fill(McParticle->Eta(), MotherParticle->Eta());
332 fHMuMotherRecEta[1]->Fill(muonTrack->Eta(), MotherParticle->Eta());
333 fHMuDCA[1]->Fill(muonTrack->GetDCA());
334 }
335 else
336 {
337 fHMuMotherGenPt[3]->Fill(McParticle->Pt(), MotherParticle->Pt());
338 fHMuMotherRecPt[3]->Fill(muonTrack->Pt(), MotherParticle->Pt());
339 fHMuMotherGenPhi[3]->Fill(McParticle->Phi(), MotherParticle->Phi());
340 fHMuMotherRecPhi[3]->Fill(muonTrack->Phi(), MotherParticle->Phi());
341 fHMuMotherGenEta[3]->Fill(McParticle->Eta(), MotherParticle->Eta());
342 fHMuMotherRecEta[3]->Fill(muonTrack->Eta(), MotherParticle->Eta());
343 fHMuDCA[3]->Fill(muonTrack->GetDCA());
344 }
345 }
346 } // (mother hadron) : (daughter muon) QA
347 }
348 }
349
350 for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
351 {
352 AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(ipart);
353 if(TMath::Abs(McParticle->PdgCode())==13)
354 {
355 Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
356 fHMuonParGen->Fill(fillArrayParGen);
357 }
358 }
359 }
360 PostData(1, fOutputList);
361 return;
362}
363
364//________________________________________________________________________
365void AliMuonEffMC::Terminate(Option_t *)
366{
367 // Draw result to the screen
368 // Called once at the end of the query
369
370 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
371 if (!fOutputList) {
372 AliError("Output list not available");
373 return;
374 }
375}
376
377//________________________________________________________________________
378Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
379{
380 // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
381
382 Int_t nContributors = 0;
383 Double_t zVertex = 999;
384 TString name("");
385
386 if (obj->InheritsFrom("AliESDEvent")) {
387 AliESDEvent* esdevt = (AliESDEvent*) obj;
388 const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
389 if (!vtx)
390 return 0;
391 nContributors = vtx->GetNContributors();
392 zVertex = vtx->GetZ();
393 name = vtx->GetName();
394 }
395 else if (obj->InheritsFrom("AliAODEvent")) {
396 AliAODEvent* aodevt = (AliAODEvent*) obj;
397 if (aodevt->GetNumberOfVertices() < 1)
398 return 0;
399 const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
400 nContributors = vtx->GetNContributors();
401 zVertex = vtx->GetZ();
402 name = vtx->GetName();
403 }
404
405 // Reject if TPC-only vertex
406 if (name.CompareTo("TPCVertex")==0)
407 return kFALSE;
408
409 // Check # contributors and range...
410 if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
411 return kFALSE;
412 }
413
414 return kTRUE;
415}
416
417
418//________________________________________________________________________
419Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
420{
421 // Applying track cuts for MUON tracks
422 if(!track.ContainTrackerData()) return kFALSE;
423 if(!track.ContainTriggerData()) return kFALSE;
424
425 Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
426 Double_t eta = track.Eta();
427
428 // Theta cut at absorber end
429 if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
430 // Eta cut
431 if(eta <= -4. || eta >= -2.5) return kFALSE;
432 if(track.GetMatchTrigger() <= 0) return kFALSE;
433
434 return kTRUE;
435}