30ee9077eeb0b4de16a3ee13cf9a0b133e8f0161
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / FourierDecomposition / AliMuonEffMC.cxx
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 "TParticle.h"
15 #include "AliStack.h"
16 #include "AliAnalysisManager.h"
17 #include "AliAnalysisTask.h"
18 #include "AliESDVertex.h"
19 #include "AliESDEvent.h"
20 #include "AliAODEvent.h"
21 #include "AliESDMuonTrack.h"
22 #include "AliAODTrack.h"
23 #include "AliESDVertex.h"
24 #include "AliAODVertex.h"
25 #include "AliCentrality.h"
26 #include "AliVParticle.h"
27 #include "AliMCParticle.h"
28 #include "AliAnalysisHelperJetTasks.h"
29 #include "AliMCEvent.h"
30 #include "AliAODMCParticle.h"
31
32 using std::cout;
33 using std::endl;
34
35 ClassImp(AliMuonEffMC)
36
37 //________________________________________________________________________
38 AliMuonEffMC::AliMuonEffMC() :
39   AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
40   fHEventStat(0), fHEvt(0x0), fIsMc(kFALSE), fMDProcess(kFALSE), fCentralityEstimator("V0M"),
41   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(36),
42   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0), fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0), fHFXmuonP(0), fHFXmuonM(0)
43 {
44   // Constructor
45   //DefineInput(0, TChain::Class());
46   //DefineOutput(1, TList::Class());
47
48
49 }
50
51 //________________________________________________________________________
52 AliMuonEffMC::AliMuonEffMC(const char *name) :
53   AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),      
54   fHEventStat(0), fHEvt(0x0),  fIsMc(kFALSE), fMDProcess(kFALSE), fCentralityEstimator("V0M"),
55   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(36),
56   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0), fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0), fHFXmuonP(0), fHFXmuonM(0)
57 {
58   // Constructor
59   DefineInput(0, TChain::Class());
60   DefineOutput(1, TList::Class());
61
62 }
63
64 //________________________________________________________________________
65 AliMuonEffMC::~AliMuonEffMC()
66 {
67   //Destructor
68   if(fOutputList) delete fOutputList;
69 }
70
71 //________________________________________________________________________
72 void AliMuonEffMC::UserCreateOutputObjects()
73 {
74   // Create histograms
75   // Called once (per slave on PROOF!)
76
77   fOutputList = new TList();
78   fOutputList->SetOwner(1);
79  
80   for(Int_t i=0; i<4; i++)
81   {
82     fHMuMotherGenPt[i] = 0x0;
83     fHMuMotherRecPt[i] = 0x0;
84     fHMuMotherGenPhi[i] = 0x0;
85     fHMuMotherRecPhi[i] = 0x0;
86     fHMuMotherGenEta[i] = 0x0;
87     fHMuMotherRecEta[i] = 0x0;
88     fHMuDCA[i] = 0x0;
89     for(Int_t j=0; j<3; j++)
90     {
91       fHMuMohterPhiDifGen[i][j] = 0x0;
92       fHMuMohterPhiDifRec[i][j] = 0x0;
93       fHMuMohterEtaDifGen[i][j] = 0x0;
94       fHMuMohterEtaDifRec[i][j] = 0x0;
95     }
96   }
97
98   fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
99   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
100   fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
101   fHEventStat->GetXaxis()->SetBinLabel(3,"File");
102   fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
103   fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
104   fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
105   fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
106   fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
107   fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
108   fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
109   fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
110   fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
111   fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
112   fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
113   fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
114   fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
115   fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
116   fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
117   fOutputList->Add(fHEventStat);
118
119   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
120   fOutputList->Add(fHEvt);
121
122
123   Int_t iTrackBin[5];
124   Double_t* trackBins[5];
125   const char* trackAxisTitle[5];
126
127   // eta
128   Double_t etaBins[fNEtaBins+1];
129   for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
130   iTrackBin[0] = fNEtaBins;
131   trackBins[0] = etaBins;
132   trackAxisTitle[0] = "#eta";
133
134   // p_T
135   Double_t pTBins[fNpTBins+1];
136   for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
137   iTrackBin[1] = fNpTBins;
138   trackBins[1] = pTBins;
139   trackAxisTitle[1] = "p_{T} (GeV/c)";
140
141   // centrality
142   Double_t CentBins[fNCentBins+1];
143   for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
144   iTrackBin[2] = fNCentBins;
145   trackBins[2] = CentBins;
146   trackAxisTitle[2] = "Cent";
147
148   // Z-vertex
149   Double_t ZvtxBins[fNZvtxBins+1];
150   for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
151   iTrackBin[3] = fNZvtxBins;
152   trackBins[3] = ZvtxBins;
153   trackAxisTitle[3] = "Zvtx";
154
155   // phi
156   Double_t phiBins[fNPhiBins+1];
157   for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
158   iTrackBin[4] = fNPhiBins;
159   trackBins[4] = phiBins;
160   trackAxisTitle[4] = "#phi";
161
162   // THn for tracking efficiency
163   fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
164   for (Int_t i=0; i<5; i++)
165   {
166     fHMuonParGen->SetBinEdges(i, trackBins[i]);
167     fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
168   }
169   fHMuonParGen->Sumw2();
170   fOutputList->Add(fHMuonParGen);
171
172   fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
173   fHMuonDetGen->Sumw2();
174   fOutputList->Add(fHMuonDetGen);
175
176   fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
177   fHMuonDetRec->Sumw2();
178   fOutputList->Add(fHMuonDetRec);
179
180   fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
181   fHEtcDetRec->Sumw2();
182   fOutputList->Add(fHEtcDetRec);
183
184   if(fIsMc)
185   {
186     const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
187     const char *MuPt[3] = {"0510","1020","2040"};
188     if(fMDProcess)
189     {
190       for(Int_t i=0; i<4; i++)
191       {
192         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);
193         fOutputList->Add(fHMuMotherGenPt[i]);
194         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);
195         fOutputList->Add(fHMuMotherRecPt[i]);
196         fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
197         fOutputList->Add(fHMuMotherGenPhi[i]);
198         fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
199         fOutputList->Add(fHMuMotherRecPhi[i]);
200         fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
201         fOutputList->Add(fHMuMotherGenEta[i]);
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         fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
205         fOutputList->Add(fHMuDCA[i]);
206         
207         for(Int_t j=0; j<3; j++)
208         {
209           fHMuMohterPhiDifGen[i][j] = new TH1F(Form("fHMuMohterPhiDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
210           fOutputList->Add(fHMuMohterPhiDifGen[i][j]);
211           fHMuMohterPhiDifRec[i][j] = new TH1F(Form("fHMuMohterPhiDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
212           fOutputList->Add(fHMuMohterPhiDifRec[i][j]);
213           fHMuMohterEtaDifGen[i][j] = new TH1F(Form("fHMuMohterEtaDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
214           fOutputList->Add(fHMuMohterEtaDifGen[i][j]);
215           fHMuMohterEtaDifRec[i][j] = new TH1F(Form("fHMuMohterEtaDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
216           fOutputList->Add(fHMuMohterEtaDifRec[i][j]);
217         }
218       }
219     }
220     
221     fHFXu = new TH1F("fHFXu",";x_{F}",1000, 0.0, 1.0);
222     fOutputList->Add(fHFXu);
223     
224     fHFXantiu = new TH1F("fHFXantiu",";x_{F}",1000, 0.0, 1.0);
225     fOutputList->Add(fHFXantiu);
226     
227     fHFXd = new TH1F("fHFXd",";x_{F}",1000, 0.0, 1.0);
228     fOutputList->Add(fHFXd);
229     
230     fHFXantid = new TH1F("fHFXantid",";x_{F}",1000, 0.0, 1.0);
231     fOutputList->Add(fHFXantid);
232     
233     fHFXg = new TH1F("fHFXg",";x_{F}",1000, 0.0, 1.0);
234     fOutputList->Add(fHFXg);
235     
236     fHFXetc = new TH1F("fHFXetc",";x_{F}",1000, 0.0, 1.0);
237     fOutputList->Add(fHFXetc);
238     
239     fHFXmuonP = new TH1F("fHFXmuonP",";x_{F}",1000, 0.0, 1.0);
240     fOutputList->Add(fHFXmuonP);
241     
242     fHFXmuonM = new TH1F("fHFXmuonM",";x_{F}",1000, 0.0, 1.0);
243     fOutputList->Add(fHFXmuonM);
244   }
245   PostData(1, fOutputList);
246 }
247
248 //________________________________________________________________________
249 void AliMuonEffMC::UserExec(Option_t *) 
250 {
251   // Main loop, Called for each event
252   Int_t ntrks = 0; // number of tracks in an event
253
254   if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
255   {
256     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
257     if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
258     ntrks = fAOD->GetNTracks();
259   } 
260   else
261   {
262     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
263     if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
264     ntrks = fESD->GetNumberOfMuonTracks();
265   }
266         
267   fHEventStat->Fill(0.5);
268   if(fIsMc) 
269   {
270     fMC = MCEvent();
271   }
272
273   // Centrality, vertex, other event variables...
274   if(fAOD)
275   { 
276     const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
277     fZVertex = vertex->GetZ();
278     if(fAOD->GetCentrality())  fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
279   }
280   else if(fESD) 
281   {
282     const AliESDVertex* vertex = fESD->GetPrimaryVertex();
283     fZVertex = vertex->GetZ();
284     if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
285   }   
286
287   if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex)); 
288     return; }
289   if (fCentrality > 100. || fCentrality < -1.5) { //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); 
290     return; }
291  
292   if(fCentrality < 0) fCentrality = 0.0; //ad hoc centrality for pp
293   // Fill Event histogram
294   fHEvt->Fill(fZVertex, fCentrality);
295   fHEventStat->Fill(1.5);
296  
297   ULong64_t trigword = 0;
298   if(fAOD) trigword=fAOD->GetTriggerMask();
299   else if(fESD) trigword=fESD->GetTriggerMask();
300  
301   if (trigword & 0x01) fHEventStat->Fill(17.5);
302   if (trigword & 0x02) fHEventStat->Fill(3.5);
303   if (trigword & 0x04) fHEventStat->Fill(4.5);
304   if (trigword & 0x08) fHEventStat->Fill(5.5);      
305   if (trigword & 0x010) fHEventStat->Fill(6.5);
306   if (trigword & 0x020) fHEventStat->Fill(7.5);
307   if (trigword & 0x040) fHEventStat->Fill(8.5);
308   if (trigword & 0x080) fHEventStat->Fill(9.5);
309   if (trigword & 0x100) fHEventStat->Fill(10.5);
310   if (trigword & 0x200) fHEventStat->Fill(11.5);
311   if (trigword & 0x400) fHEventStat->Fill(12.5);
312   if (trigword & 0x800) fHEventStat->Fill(13.5);
313   if (trigword & 0x1000) fHEventStat->Fill(14.5);
314   if (trigword & 0x2000) fHEventStat->Fill(15.5);
315   if (trigword & 0x4000) fHEventStat->Fill(16.5);
316   
317   // generated level loop
318   if(fIsMc)
319   {
320     for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
321     {
322       if(fAOD)
323       {
324         AliAODMCParticle *AodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(ipart);
325         if(TMath::Abs(AodMcParticle->PdgCode())==13) 
326         {
327           Double_t fillArrayParGen[5] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi() };
328           fHMuonParGen->Fill(fillArrayParGen);
329         }
330       } 
331       else if(fESD)
332       {
333         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
334         if(TMath::Abs(McParticle->PdgCode())==13) 
335         {
336           Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
337           fHMuonParGen->Fill(fillArrayParGen);
338         }
339       }
340     }
341   }
342   
343   // reconstructed level loop
344   for (Int_t iTrack = 0; iTrack<ntrks; iTrack++) 
345   { 
346     Int_t label = 0;
347     Double_t trackpt = 0;
348     Double_t tracketa = 0;
349     Double_t trackphi = 0;
350     Double_t dcavalue = 0;
351     
352     Double_t mcpt = 0;
353     Double_t mceta = 0;
354     Double_t mcphi = 0;
355     
356     Int_t motherlabel =0;
357     Int_t motherpdg = 0;
358     Double_t motherpt = 0;
359     Double_t mothereta = 0;
360     Double_t motherphi = 0;
361     
362     if(fAOD)
363     {
364       AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
365       if(muonTrack)
366       {
367         if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
368         trackpt = muonTrack->Pt();
369         tracketa = muonTrack->Eta();
370         trackphi = muonTrack->Phi();
371         label =  TMath::Abs(muonTrack->GetLabel());
372       }
373     }
374     else if(fESD)
375     {
376       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
377       if(muonTrack)
378       {
379         if(!IsGoodMUONtrack(*muonTrack)) continue;
380         trackpt = muonTrack->Pt();
381         tracketa = muonTrack->Eta();
382         trackphi = muonTrack->Phi();
383         label =  TMath::Abs(muonTrack->GetLabel());
384         dcavalue = muonTrack->GetDCA();
385       }
386     }
387     Double_t fillArrayDetRec[5] = { tracketa, trackpt, fCentrality, fZVertex, trackphi };
388     fHMuonDetRec->Fill(fillArrayDetRec); 
389     
390     if(fIsMc)
391     {
392       if(fAOD)
393       {
394         AliAODMCParticle *aodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(label);
395         if(TMath::Abs(aodMcParticle->PdgCode()) != 13) 
396         {
397           fHEtcDetRec->Fill(fillArrayDetRec); 
398           continue;
399         }
400         mcpt = aodMcParticle->Pt();
401         mceta = aodMcParticle->Eta();
402         mcphi = aodMcParticle->Phi();
403         motherlabel = aodMcParticle->GetMother();
404         if(motherlabel > 0) 
405         { 
406           AliAODMCParticle *aodMotherParticle  = (AliAODMCParticle*)fMC->GetTrack(motherlabel);
407           motherpdg = TMath::Abs(aodMotherParticle->PdgCode());
408           motherpt = aodMotherParticle->Pt();
409           mothereta = aodMotherParticle->Eta();
410           motherphi = aodMotherParticle->Phi();
411         }      
412       }
413       else if(fESD)
414       {
415         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
416         if(TMath::Abs(McParticle->PdgCode()) != 13) 
417         {
418           fHEtcDetRec->Fill(fillArrayDetRec); 
419           continue;
420         }
421         mcpt = McParticle->Pt();
422         mceta = McParticle->Eta();
423         mcphi = McParticle->Phi();
424         motherlabel = McParticle->GetMother();
425         if(motherlabel > 0) 
426         { 
427           AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
428           motherpdg = TMath::Abs(MotherParticle->PdgCode());
429           motherpt = MotherParticle->Pt();
430           mothereta = MotherParticle->Eta();
431           motherphi = MotherParticle->Phi();
432         }
433       }
434       Double_t fillArrayDetGen[5] = { mceta, mcpt, fCentrality, fZVertex, mcphi };
435       fHMuonDetGen->Fill(fillArrayDetGen);
436     
437       // mother-daughter kinematic relation
438       if(fMDProcess && motherlabel>0) MDProcess(motherpdg, mcpt, mcphi, mceta, trackpt, trackphi, tracketa, motherpt, motherphi, mothereta, dcavalue);
439     }
440   }
441   PostData(1, fOutputList);
442   return;
443 }
444
445 //________________________________________________________________________
446 void AliMuonEffMC::Terminate(Option_t *) 
447 {
448   // Draw result to the screen
449   // Called once at the end of the query
450 }
451 //________________________________________________________________________
452 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
453 {
454   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
455   
456   Int_t nContributors  = 0;
457   Double_t zVertex     = 999;
458   TString name("");
459   
460   if (obj->InheritsFrom("AliESDEvent")) {
461     AliESDEvent* esdevt = (AliESDEvent*) obj;
462     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
463     if (!vtx)
464       return 0;
465     nContributors = vtx->GetNContributors();
466     zVertex       = vtx->GetZ();
467     name          = vtx->GetName();
468   }
469   else if (obj->InheritsFrom("AliAODEvent")) { 
470     AliAODEvent* aodevt = (AliAODEvent*) obj;
471     if (aodevt->GetNumberOfVertices() < 1)
472       return 0;
473     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
474     nContributors = vtx->GetNContributors();
475     zVertex       = vtx->GetZ();
476     name          = vtx->GetName();
477   }
478   
479   // Reject if TPC-only vertex
480   if (name.CompareTo("TPCVertex")==0)
481     return kFALSE;
482   
483   // Check # contributors and range...
484   if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
485     return kFALSE;
486   }
487   
488   return kTRUE;
489 }
490
491
492 //________________________________________________________________________
493 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
494 {
495   // Applying track cuts for MUON tracks
496   if(!track.ContainTrackerData()) return kFALSE;
497
498   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
499   Double_t eta = track.Eta();
500
501   // Theta cut at absorber end
502   if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
503   // Eta cut
504   if(eta <= -4. || eta >= -2.5) return kFALSE;
505   if(track.GetMatchTrigger() <= 0) return kFALSE;
506
507   return kTRUE;
508 }
509
510 //________________________________________________________________________
511 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
512 {
513   if (!track.IsMuonTrack()) return kFALSE;
514
515   Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)
516                      * TMath::RadToDeg();
517   if ((dThetaAbs<2.) || (dThetaAbs>10.)) return kFALSE;
518
519   Double_t dEta = track.Eta();
520   if ((dEta<-4.) || (dEta>2.5)) return kFALSE;
521
522   if (track.GetMatchTrigger()<0.5) return kFALSE;
523
524   return kTRUE;
525 }
526
527 //________________________________________________________________________
528 void AliMuonEffMC::MDProcess(Int_t motherpdg, Double_t mcpt, Double_t mcphi, Double_t mceta, Double_t trackpt, Double_t trackphi, Double_t tracketa, Double_t motherpt, Double_t motherphi, Double_t mothereta, Double_t dcavalue)
529 {
530   if(motherpdg==411 || motherpdg==413 || motherpdg==421 || motherpdg==423 || motherpdg==431 || motherpdg==433 || motherpdg==10413 || motherpdg==10411 || motherpdg==10423 || motherpdg==10421 || motherpdg==10433 || motherpdg==10431 || motherpdg==20413 || motherpdg==415 || motherpdg==20423 || motherpdg==425 || motherpdg==20433 || motherpdg==435)
531   {  
532     fHMuMotherGenPt[2]->Fill(mcpt, motherpt);
533     fHMuMotherRecPt[2]->Fill(trackpt, motherpt);
534     fHMuMotherGenPhi[2]->Fill(mcphi, motherphi);
535     fHMuMotherRecPhi[2]->Fill(trackphi, motherphi);
536     fHMuMotherGenEta[2]->Fill(mceta, mothereta);
537     fHMuMotherRecEta[2]->Fill(tracketa, mothereta);
538     
539     // generated 
540     if((0.5 <= mcpt) && (mcpt < 1.0))
541     {
542       fHMuMohterPhiDifGen[2][0]->Fill(deltaphi(motherphi-mcphi));
543       fHMuMohterEtaDifGen[2][0]->Fill(mothereta-mceta);
544     }
545     else if((1.0 <= mcpt) && (mcpt < 2.0))
546     {
547       fHMuMohterPhiDifGen[2][1]->Fill(deltaphi(motherphi-mcphi));
548       fHMuMohterEtaDifGen[2][1]->Fill(mothereta-mceta);
549     }
550     else if((2.0 <= mcpt) && (mcpt < 4.0))
551     {
552       fHMuMohterPhiDifGen[2][2]->Fill(deltaphi(motherphi-mcphi));
553       fHMuMohterEtaDifGen[2][2]->Fill(mothereta-mceta);
554     }
555     // reconstructed
556     if((0.5 <= trackpt) && (trackpt < 1.0))
557     {
558       fHMuMohterPhiDifRec[2][0]->Fill(deltaphi(motherphi-trackphi));
559       fHMuMohterEtaDifRec[2][0]->Fill(mothereta-tracketa);
560     }
561     else if((1.0 <= trackpt) && (trackpt < 2.0))
562     {
563       fHMuMohterPhiDifRec[2][1]->Fill(deltaphi(motherphi-trackphi));
564       fHMuMohterEtaDifRec[2][1]->Fill(mothereta-tracketa);
565     }
566     else if((2.0 <= trackpt) && (trackpt < 4.0))
567     {
568       fHMuMohterPhiDifRec[2][2]->Fill(deltaphi(motherphi-trackphi));
569       fHMuMohterEtaDifRec[2][2]->Fill(mothereta-tracketa);
570     }
571     
572     // DCA
573     if(fESD) fHMuDCA[2]->Fill(dcavalue);
574   }
575   
576   else if(motherpdg==211) 
577   {
578     fHMuMotherGenPt[0]->Fill(mcpt, motherpt);
579     fHMuMotherRecPt[0]->Fill(trackpt, motherpt);
580     fHMuMotherGenPhi[0]->Fill(mcphi, motherphi);
581     fHMuMotherRecPhi[0]->Fill(trackphi, motherphi);
582     fHMuMotherGenEta[0]->Fill(mceta, mothereta);
583     fHMuMotherRecEta[0]->Fill(tracketa, mothereta);
584     // generated 
585     if((0.5 <= mcpt) && (mcpt < 1.0))
586     {
587       fHMuMohterPhiDifGen[0][0]->Fill(deltaphi(motherphi-mcphi));
588       fHMuMohterEtaDifGen[0][0]->Fill(mothereta-mceta);
589     }
590     else if((1.0 <= mcpt) && (mcpt < 2.0))
591     {
592       fHMuMohterPhiDifGen[0][1]->Fill(deltaphi(motherphi-mcphi));
593       fHMuMohterEtaDifGen[0][1]->Fill(mothereta-mceta);
594     }
595     else if((2.0 <= mcpt) && (mcpt < 4.0))
596     {
597       fHMuMohterPhiDifGen[0][2]->Fill(deltaphi(motherphi-mcphi));
598       fHMuMohterEtaDifGen[0][2]->Fill(mothereta-mceta);
599     }
600     // reconstructed
601     if((0.5 <= trackpt) && (trackpt < 1.0))
602     {
603       fHMuMohterPhiDifRec[0][0]->Fill(deltaphi(motherphi-trackphi));
604       fHMuMohterEtaDifRec[0][0]->Fill(mothereta-tracketa);
605     }
606     else if((1.0 <= trackpt) && (trackpt < 2.0))
607     {
608       fHMuMohterPhiDifRec[0][1]->Fill(deltaphi(motherphi-trackphi));
609       fHMuMohterEtaDifRec[0][1]->Fill(mothereta-tracketa);
610     }
611     else if((2.0 <= trackpt) && (trackpt < 4.0))
612     {
613       fHMuMohterPhiDifRec[0][2]->Fill(deltaphi(motherphi-trackphi));
614       fHMuMohterEtaDifRec[0][2]->Fill(mothereta-tracketa);
615     }
616     if(fESD) fHMuDCA[0]->Fill(dcavalue);
617   }
618   else if(motherpdg==321)
619   {
620     fHMuMotherGenPt[1]->Fill(mcpt, motherpt);
621     fHMuMotherRecPt[1]->Fill(trackpt, motherpt);
622     fHMuMotherGenPhi[1]->Fill(mcphi, motherphi);
623     fHMuMotherRecPhi[1]->Fill(trackphi, motherphi);
624     fHMuMotherGenEta[1]->Fill(mceta, mothereta);
625     fHMuMotherRecEta[1]->Fill(tracketa, mothereta);
626     // generated 
627     if((0.5 <= mcpt) && (mcpt < 1.0))
628     {
629       fHMuMohterPhiDifGen[1][0]->Fill(deltaphi(motherphi-mcphi));
630       fHMuMohterEtaDifGen[1][0]->Fill(mothereta-mceta);
631     }
632     else if((1.0 <= mcpt) && (mcpt < 2.0))
633     {
634       fHMuMohterPhiDifGen[1][1]->Fill(deltaphi(motherphi-mcphi));
635       fHMuMohterEtaDifGen[1][1]->Fill(mothereta-mceta);
636     }
637     else if((2.0 <= mcpt) && (mcpt < 4.0))
638     {
639       fHMuMohterPhiDifGen[1][2]->Fill(deltaphi(motherphi-mcphi));
640       fHMuMohterEtaDifGen[1][2]->Fill(mothereta-mceta);
641     }
642     // reconstructed
643     if((0.5 <= trackpt) && (trackpt < 1.0))
644     {
645       fHMuMohterPhiDifRec[1][0]->Fill(deltaphi(motherphi-trackphi));
646       fHMuMohterEtaDifRec[1][0]->Fill(mothereta-tracketa);
647     }
648     else if((1.0 <= trackpt) && (trackpt < 2.0))
649     {
650       fHMuMohterPhiDifRec[1][1]->Fill(deltaphi(motherphi-trackphi));
651       fHMuMohterEtaDifRec[1][1]->Fill(mothereta-tracketa);
652     }
653     else if((2.0 <= trackpt) && (trackpt < 4.0))
654     {
655       fHMuMohterPhiDifRec[1][2]->Fill(deltaphi(motherphi-trackphi));
656       fHMuMohterEtaDifRec[1][2]->Fill(mothereta-tracketa);
657     }
658     if(fESD) fHMuDCA[1]->Fill(dcavalue);
659   }     
660   else 
661   {
662     fHMuMotherGenPt[3]->Fill(mcpt, motherpt);
663     fHMuMotherRecPt[3]->Fill(trackpt, motherpt);
664     fHMuMotherGenPhi[3]->Fill(mcphi, motherphi);
665     fHMuMotherRecPhi[3]->Fill(trackphi, motherphi);
666     fHMuMotherGenEta[3]->Fill(mceta, mothereta);
667     fHMuMotherRecEta[3]->Fill(tracketa, mothereta);
668     // generated 
669     if((0.5 <= mcpt) && (mcpt < 1.0))
670     {
671       fHMuMohterPhiDifGen[3][0]->Fill(deltaphi(motherphi-mcphi));
672       fHMuMohterEtaDifGen[3][0]->Fill(mothereta-mceta);
673     }
674     else if((1.0 <= mcpt) && (mcpt < 2.0))
675     {
676       fHMuMohterPhiDifGen[3][1]->Fill(deltaphi(motherphi-mcphi));
677       fHMuMohterEtaDifGen[3][1]->Fill(mothereta-mceta);
678     }
679     else if((2.0 <= mcpt) && (mcpt < 4.0))
680     {
681       fHMuMohterPhiDifGen[3][2]->Fill(deltaphi(motherphi-mcphi));
682       fHMuMohterEtaDifGen[3][2]->Fill(mothereta-mceta);
683     }
684     // reconstructed
685     if((0.5 <= trackpt) && (trackpt < 1.0))
686     {
687       fHMuMohterPhiDifRec[3][0]->Fill(deltaphi(motherphi-trackphi));
688       fHMuMohterEtaDifRec[3][0]->Fill(mothereta-tracketa);
689     }
690     else if((1.0 <= trackpt) && (trackpt < 2.0))
691     {
692       fHMuMohterPhiDifRec[3][1]->Fill(deltaphi(motherphi-trackphi));
693       fHMuMohterEtaDifRec[3][1]->Fill(mothereta-tracketa);
694     }
695     else if((2.0 <= trackpt) && (trackpt < 4.0))
696     {
697       fHMuMohterPhiDifRec[3][2]->Fill(deltaphi(motherphi-trackphi));
698       fHMuMohterEtaDifRec[3][2]->Fill(mothereta-tracketa);
699     }
700     if(fESD) fHMuDCA[3]->Fill(dcavalue);
701   }  
702 }
703 //________________________________________________________________________
704 Double_t AliMuonEffMC::deltaphi(Double_t phi)
705 {
706   if(phi < -1.0*TMath::Pi()) { return (phi + TMath::TwoPi()); }
707   else if(phi > TMath::Pi()) { return (phi - TMath::TwoPi()); }
708   else { return phi; }
709 }