f5b04a6145501509effc494ab1893326f5d62f33
[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 "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
28 using std::cout;
29 using std::endl;
30
31 ClassImp(AliMuonEffMC)
32
33 //________________________________________________________________________
34 AliMuonEffMC::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 //________________________________________________________________________
57 AliMuonEffMC::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 //________________________________________________________________________
80 AliMuonEffMC::~AliMuonEffMC()
81 {
82   //Destructor
83 }
84
85 //________________________________________________________________________
86 void 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 //________________________________________________________________________
212 void 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 //________________________________________________________________________
365 void 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 //________________________________________________________________________
378 Bool_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 //________________________________________________________________________
419 Bool_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 }