1 // MUON track QA referring AliMuonEffMC.cxx
2 // Author : Saehanseul Oh
4 #include "AliMuonEffMC.h"
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 "AliAODTrack.h"
21 #include "AliESDVertex.h"
22 #include "AliAODVertex.h"
23 #include "AliCentrality.h"
24 #include "AliVParticle.h"
25 #include "AliMCParticle.h"
26 #include "AliAnalysisHelperJetTasks.h"
27 #include "AliMCEvent.h"
32 ClassImp(AliMuonEffMC)
34 //________________________________________________________________________
35 AliMuonEffMC::AliMuonEffMC() :
36 AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
37 fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
38 fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
39 fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
44 //________________________________________________________________________
45 AliMuonEffMC::AliMuonEffMC(const char *name) :
46 AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
47 fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
48 fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
49 fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
52 DefineInput(0, TChain::Class());
53 DefineOutput(1, TList::Class());
56 //________________________________________________________________________
57 AliMuonEffMC::~AliMuonEffMC()
60 if(fOutputList) delete fOutputList;
63 //________________________________________________________________________
64 void AliMuonEffMC::UserCreateOutputObjects()
67 // Called once (per slave on PROOF!)
68 for(Int_t i=0; i<4; i++)
70 fHMuMotherGenPt[i] = 0x0;
71 fHMuMotherRecPt[i] = 0x0;
72 fHMuMotherGenPhi[i] = 0x0;
73 fHMuMotherRecPhi[i] = 0x0;
74 fHMuMotherGenEta[i] = 0x0;
75 fHMuMotherRecEta[i] = 0x0;
79 fOutputList = new TList();
80 fOutputList->SetOwner(1);
82 fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
83 fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
84 fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
85 fHEventStat->GetXaxis()->SetBinLabel(3,"File");
86 fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt"); //!Global Trigger Single plus High p_T
87 fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt"); //!Global Trigger Single plus All p_T
88 fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt"); //!Global Trigger Single minus Low p_T
89 fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt"); //!Global Trigger Single minus High p_T
90 fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt"); //!Global Trigger Single minus All p_T
91 fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt"); //!Global Trigger Single undefined Low p_T
92 fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
93 fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt"); //!Global Trigger Single undefined All p_T
94 fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt"); //!Global Trigger UnlikeSign pair Low p_T
95 fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
96 fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt"); //!Global Trigger UnlikeSign pair All p_T
97 fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt"); //!Global Trigger LikeSign pair pair Low p_T
98 fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
99 fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt"); //!Global Trigger LikeSign pair pair All p_T
100 fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt"); //!Global Trigger Single plus Low p_T
101 fOutputList->Add(fHEventStat);
103 fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
104 fOutputList->Add(fHEvt);
108 Double_t* trackBins[5];
109 const char* trackAxisTitle[5];
112 Double_t etaBins[fNEtaBins+1];
113 for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
114 iTrackBin[0] = fNEtaBins;
115 trackBins[0] = etaBins;
116 trackAxisTitle[0] = "#eta";
119 Double_t pTBins[fNpTBins+1];
120 for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
121 iTrackBin[1] = fNpTBins;
122 trackBins[1] = pTBins;
123 trackAxisTitle[1] = "p_{T} (GeV/c)";
126 Double_t CentBins[fNCentBins+1];
127 for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
128 iTrackBin[2] = fNCentBins;
129 trackBins[2] = CentBins;
130 trackAxisTitle[2] = "Cent";
133 Double_t ZvtxBins[fNZvtxBins+1];
134 for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
135 iTrackBin[3] = fNZvtxBins;
136 trackBins[3] = ZvtxBins;
137 trackAxisTitle[3] = "Zvtx";
140 Double_t phiBins[fNPhiBins+1];
141 for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
142 iTrackBin[4] = fNPhiBins;
143 trackBins[4] = phiBins;
144 trackAxisTitle[4] = "#phi";
146 // THn for tracking efficiency
147 fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
148 for (Int_t i=0; i<5; i++)
150 fHMuonParGen->SetBinEdges(i, trackBins[i]);
151 fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
153 fHMuonParGen->Sumw2();
154 fOutputList->Add(fHMuonParGen);
156 fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
157 fHMuonDetGen->Sumw2();
158 fOutputList->Add(fHMuonDetGen);
160 fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
161 fHMuonDetRec->Sumw2();
162 fOutputList->Add(fHMuonDetRec);
164 fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
165 fHEtcDetRec->Sumw2();
166 fOutputList->Add(fHEtcDetRec);
170 const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
172 for(Int_t i=0; i<4; i++)
174 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);
175 fOutputList->Add(fHMuMotherGenPt[i]);
177 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);
178 fOutputList->Add(fHMuMotherRecPt[i]);
180 fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
181 fOutputList->Add(fHMuMotherGenPhi[i]);
183 fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
184 fOutputList->Add(fHMuMotherRecPhi[i]);
186 fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
187 fOutputList->Add(fHMuMotherGenEta[i]);
189 fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
190 fOutputList->Add(fHMuMotherRecEta[i]);
192 fHMuDCA[i] = new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
193 fOutputList->Add(fHMuDCA[i]);
196 PostData(1, fOutputList);
199 //________________________________________________________________________
200 void AliMuonEffMC::UserExec(Option_t *)
202 // Main loop, Called for each event
203 Int_t ntrks = 0; // number of tracks in an event
205 if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
207 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
208 if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
209 ntrks = fAOD->GetNTracks();
213 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
214 if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
215 ntrks = fESD->GetNumberOfMuonTracks();
218 fHEventStat->Fill(0.5);
222 if (!fMC) { AliError("MC event not avaliable."); return; }
225 // Centrality, vertex, other event variables...
228 const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
229 fZVertex = vertex->GetZ();
230 if(fAOD->GetCentrality()) fCentrality = fAOD->GetCentrality()->GetCentralityPercentile("fCentralityEstimator");
234 const AliESDVertex* vertex = fESD->GetPrimaryVertex();
235 fZVertex = vertex->GetZ();
236 if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile("fCentralityEstimator");
239 if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { AliInfo(Form("Event REJECTED. z = %.1f", fZVertex)); return; }
240 if (fCentrality > 100. || fCentrality < -1.5) { AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); return; }
242 if(fCentrality < 0) fCentrality = 100.5;
243 // Fill Event histogram
244 fHEvt->Fill(fZVertex, fCentrality);
245 fHEventStat->Fill(1.5);
247 ULong64_t trigword = 0;
248 if(fAOD) trigword=fAOD->GetTriggerMask();
249 else if(fESD) trigword=fESD->GetTriggerMask();
251 if (trigword & 0x01) fHEventStat->Fill(17.5);
252 if (trigword & 0x02) fHEventStat->Fill(3.5);
253 if (trigword & 0x04) fHEventStat->Fill(4.5);
254 if (trigword & 0x08) fHEventStat->Fill(5.5);
255 if (trigword & 0x010) fHEventStat->Fill(6.5);
256 if (trigword & 0x020) fHEventStat->Fill(7.5);
257 if (trigword & 0x040) fHEventStat->Fill(8.5);
258 if (trigword & 0x080) fHEventStat->Fill(9.5);
259 if (trigword & 0x100) fHEventStat->Fill(10.5);
260 if (trigword & 0x200) fHEventStat->Fill(11.5);
261 if (trigword & 0x400) fHEventStat->Fill(12.5);
262 if (trigword & 0x800) fHEventStat->Fill(13.5);
263 if (trigword & 0x1000) fHEventStat->Fill(14.5);
264 if (trigword & 0x2000) fHEventStat->Fill(15.5);
265 if (trigword & 0x4000) fHEventStat->Fill(16.5);
269 for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
271 // loop on muon tracks
273 Double_t trackpt = 0;
274 Double_t tracketa = 0;
275 Double_t trackphi = 0;
276 Double_t dcavalue = 0;
279 AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
282 if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
283 trackpt = muonTrack->Pt();
284 tracketa = muonTrack->Eta();
285 trackphi = muonTrack->Phi();
286 label = TMath::Abs(muonTrack->GetLabel());
291 AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
294 if(!IsGoodMUONtrack(*muonTrack)) continue;
295 trackpt = muonTrack->Pt();
296 tracketa = muonTrack->Eta();
297 trackphi = muonTrack->Phi();
298 label = TMath::Abs(muonTrack->GetLabel());
299 dcavalue = muonTrack->GetDCA();
303 Double_t fillArrayDetRec[5] = { tracketa, trackpt, fCentrality, fZVertex, trackphi };
304 fHMuonDetRec->Fill(fillArrayDetRec);
306 AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(label);
307 if(TMath::Abs(McParticle->PdgCode()) != 13)
309 fHEtcDetRec->Fill(fillArrayDetRec);
312 Double_t fillArrayDetGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
313 fHMuonDetGen->Fill(fillArrayDetGen);
317 if(McParticle->GetMother() > 0)
319 AliMCParticle *MotherParticle = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
320 Int_t motherlabel = TMath::Abs(MotherParticle->PdgCode());
322 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)
324 fHMuMotherGenPt[2]->Fill(McParticle->Pt(), MotherParticle->Pt());
325 fHMuMotherRecPt[2]->Fill(trackpt, MotherParticle->Pt());
326 fHMuMotherGenPhi[2]->Fill(McParticle->Phi(), MotherParticle->Phi());
327 fHMuMotherRecPhi[2]->Fill(trackphi, MotherParticle->Phi());
328 fHMuMotherGenEta[2]->Fill(McParticle->Eta(), MotherParticle->Eta());
329 fHMuMotherRecEta[2]->Fill(tracketa, MotherParticle->Eta());
330 if(fESD) fHMuDCA[2]->Fill(dcavalue);
332 else if(motherlabel==211)
334 fHMuMotherGenPt[0]->Fill(McParticle->Pt(), MotherParticle->Pt());
335 fHMuMotherRecPt[0]->Fill(trackpt, MotherParticle->Pt());
336 fHMuMotherGenPhi[0]->Fill(McParticle->Phi(), MotherParticle->Phi());
337 fHMuMotherRecPhi[0]->Fill(trackphi, MotherParticle->Phi());
338 fHMuMotherGenEta[0]->Fill(McParticle->Eta(), MotherParticle->Eta());
339 fHMuMotherRecEta[0]->Fill(tracketa, MotherParticle->Eta());
340 if(fESD) fHMuDCA[0]->Fill(dcavalue);
342 else if(motherlabel==321)
344 fHMuMotherGenPt[1]->Fill(McParticle->Pt(), MotherParticle->Pt());
345 fHMuMotherRecPt[1]->Fill(trackpt, MotherParticle->Pt());
346 fHMuMotherGenPhi[1]->Fill(McParticle->Phi(), MotherParticle->Phi());
347 fHMuMotherRecPhi[1]->Fill(trackphi, MotherParticle->Phi());
348 fHMuMotherGenEta[1]->Fill(McParticle->Eta(), MotherParticle->Eta());
349 fHMuMotherRecEta[1]->Fill(tracketa, MotherParticle->Eta());
350 if(fESD) fHMuDCA[1]->Fill(dcavalue);
354 fHMuMotherGenPt[3]->Fill(McParticle->Pt(), MotherParticle->Pt());
355 fHMuMotherRecPt[3]->Fill(trackpt, MotherParticle->Pt());
356 fHMuMotherGenPhi[3]->Fill(McParticle->Phi(), MotherParticle->Phi());
357 fHMuMotherRecPhi[3]->Fill(trackphi, MotherParticle->Phi());
358 fHMuMotherGenEta[3]->Fill(McParticle->Eta(), MotherParticle->Eta());
359 fHMuMotherRecEta[3]->Fill(tracketa, MotherParticle->Eta());
360 if(fESD) fHMuDCA[3]->Fill(dcavalue);
363 } // (mother hadron) : (daughter muon) QA
367 for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
369 AliMCParticle *McParticle = (AliMCParticle*)fMC->GetTrack(ipart);
370 if(TMath::Abs(McParticle->PdgCode())==13)
372 Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
373 fHMuonParGen->Fill(fillArrayParGen);
377 PostData(1, fOutputList);
381 //________________________________________________________________________
382 void AliMuonEffMC::Terminate(Option_t *)
384 // Draw result to the screen
385 // Called once at the end of the query
389 //________________________________________________________________________
390 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
392 // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
394 Int_t nContributors = 0;
395 Double_t zVertex = 999;
398 if (obj->InheritsFrom("AliESDEvent")) {
399 AliESDEvent* esdevt = (AliESDEvent*) obj;
400 const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
403 nContributors = vtx->GetNContributors();
404 zVertex = vtx->GetZ();
405 name = vtx->GetName();
407 else if (obj->InheritsFrom("AliAODEvent")) {
408 AliAODEvent* aodevt = (AliAODEvent*) obj;
409 if (aodevt->GetNumberOfVertices() < 1)
411 const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
412 nContributors = vtx->GetNContributors();
413 zVertex = vtx->GetZ();
414 name = vtx->GetName();
417 // Reject if TPC-only vertex
418 if (name.CompareTo("TPCVertex")==0)
421 // Check # contributors and range...
422 if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
430 //________________________________________________________________________
431 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
433 // Applying track cuts for MUON tracks
434 if(!track.ContainTrackerData()) return kFALSE;
436 Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
437 Double_t eta = track.Eta();
439 // Theta cut at absorber end
440 if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
442 if(eta <= -4. || eta >= -2.5) return kFALSE;
443 if(track.GetMatchTrigger() <= 0) return kFALSE;
448 //________________________________________________________________________
449 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
451 if (!track.IsMuonTrack()) return kFALSE;
453 Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)
455 if ((dThetaAbs<2.) || (dThetaAbs>10.)) return kFALSE;
457 Double_t dEta = track.Eta();
458 if ((dEta<-4.) || (dEta>2.5)) return kFALSE;
460 if (track.GetMatchTrigger()<0.5) return kFALSE;