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