9636648dd69a3b43457c1204cc6fede4e751e22e
[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   if(fOutputList) delete fOutputList;
84 }
85
86 //________________________________________________________________________
87 void 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 //________________________________________________________________________
213 void 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
245   Int_t iVerb = 0;
246
247    // Centrality, vertex, other event variables...
248   if (!VertexOk(fESD)) { 
249     if (iVerb>1)
250       AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));  
251     return; 
252   }
253
254   if (fCentrality > 100. || fCentrality < -1.5) { 
255     if (iVerb>1)
256       AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality)); 
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 //________________________________________________________________________
370 void 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 //________________________________________________________________________
383 Bool_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 //________________________________________________________________________
424 Bool_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 }