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