AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskCharmBaryonsMC.cxx
Go to the documentation of this file.
1 
2 //----------------------------------------------------------------
3 // Implementation of Class AliAnalysisTaskCharmBaryonsMC
4 //
5 // Task used to analize simulations at generation level (i.e. only
6 // needs galice.root and Kinematics.root).
7 //
8 // This task make pt and rapidity histograms of Lc/D0/Xic0
9 //
10 // Author: Yosuke Watanabe
11 //
12 //----------------------------------------------------------------
13 
14 
15 #include "TChain.h"
16 #include "TTree.h"
17 #include "TH1F.h"
18 #include "TCanvas.h"
19 
20 #include "AliAnalysisTask.h"
21 #include "AliAnalysisManager.h"
22 
23 #include "AliVEvent.h"
24 #include "AliESDEvent.h"
25 #include "AliMCEvent.h"
26 #include "AliAODEvent.h"
27 
29 #include "TGraphErrors.h"
30 #include "AliLog.h"
31 #include "AliStack.h"
32 #include "AliGenEventHeader.h"
33 #include "AliGenPythiaEventHeader.h"
34 
35 #include <iostream>
36 #include "AliPDG.h"
37 #include "AliGenDPMjetEventHeader.h"
38 #include "TFile.h"
39 
40 using namespace std;
41 
43 
45 : AliAnalysisTaskSE(), fMyOut(0),
46  fHistEvt(0),
47  fHistMult(0),
48  fHistPtPromptD0(0),
49  fHistPtPromptLc(0),
50  fHistPtPromptXic0(0),
51  fHistPtFeeddownD0(0),
52  fHistPtFeeddownLc(0),
53  fHistPtFeeddownXic0(0),
54  fHistPtInclusiveD0(0),
55  fHistPtInclusiveLc(0),
56  fHistPtInclusiveXic0(0),
57  fHistPtvsRapidityPromptD0(0),
58  fHistPtvsRapidityPromptLc(0),
59  fHistPtvsRapidityPromptXic0(0),
60  fHistPtvsRapidityFeeddownD0(0),
61  fHistPtvsRapidityFeeddownLc(0),
62  fHistPtvsRapidityFeeddownXic0(0),
63  fHistPtvsRapidityInclusiveD0(0),
64  fHistPtvsRapidityInclusiveLc(0),
65  fHistPtvsRapidityInclusiveXic0(0)
66 {
67 
68  //
69  // default constructor
70  //
71 }
72 
73 //________________________________________________________________________
75 : AliAnalysisTaskSE(name), fMyOut(0),
76  fHistEvt(0),
77  fHistMult(0),
78  fHistPtPromptD0(0),
79  fHistPtPromptLc(0),
80  fHistPtPromptXic0(0),
81  fHistPtFeeddownD0(0),
82  fHistPtFeeddownLc(0),
83  fHistPtFeeddownXic0(0),
84  fHistPtInclusiveD0(0),
85  fHistPtInclusiveLc(0),
86  fHistPtInclusiveXic0(0),
87  fHistPtvsRapidityPromptD0(0),
88  fHistPtvsRapidityPromptLc(0),
89  fHistPtvsRapidityPromptXic0(0),
90  fHistPtvsRapidityFeeddownD0(0),
91  fHistPtvsRapidityFeeddownLc(0),
92  fHistPtvsRapidityFeeddownXic0(0),
93  fHistPtvsRapidityInclusiveD0(0),
94  fHistPtvsRapidityInclusiveLc(0),
95  fHistPtvsRapidityInclusiveXic0(0)
96 {
97  //
98  // constructor
99  //
100 
101 
102  AliPDG::AddParticlesToPdgDataBase();
103 
104  DefineOutput(1, TList::Class());
105 }
106 
108 
109  // destructor
110 
111  if(fMyOut) {
112  // fMyOut owns the histos
113  delete fMyOut;
114  fMyOut = 0;
115  }
116 
117 }
118 
119 
120 //________________________________________________________________________
122 {
123  // Create histograms
124  // Called once
125 
126  fMyOut = new TList();
127  BookHistograms();
128 
129  fMyOut->SetOwner();
130 
131  // Suppress annoying printout
132  AliLog::SetGlobalLogLevel(AliLog::kError);
133 
134 }
135 
137  //
138  //
139  //
140  fHistEvt = new TH1F("fHistEvt","Number of events", 1, -0.5, 0.5);
141  fMyOut->Add(fHistEvt);
142  fHistMult = new TH1F("fHistMult","Multiplicity", 100, 0.0, 1000.);
143  fMyOut->Add(fHistMult);
144 
145  fHistPtPromptD0 = new TH1D("fHistPtPromptD0","Prompt D0", 100, 0., 20.);
146  fMyOut->Add(fHistPtPromptD0);
147  fHistPtPromptLc = new TH1D("fHistPtPromptLc","Prompt Lc", 100, 0., 20.);
148  fMyOut->Add(fHistPtPromptLc);
149  fHistPtPromptXic0 = new TH1D("fHistPtPromptXic0","Prompt Xic0", 100, 0., 20.);
151  fHistPtFeeddownD0 = new TH1D("fHistPtFeeddownD0","Feeddown D0", 100, 0., 20.);
153  fHistPtFeeddownLc = new TH1D("fHistPtFeeddownLc","Feeddown Lc", 100, 0., 20.);
155  fHistPtFeeddownXic0 = new TH1D("fHistPtFeeddownXic0","Feeddown Xic0", 100, 0., 20.);
157  fHistPtInclusiveD0 = new TH1D("fHistPtInclusiveD0","Inclusive D0", 100, 0., 20.);
159  fHistPtInclusiveLc = new TH1D("fHistPtInclusiveLc","Inclusive Lc", 100, 0., 20.);
161  fHistPtInclusiveXic0 = new TH1D("fHistPtInclusiveXic0","Inclusive Xic0", 100, 0., 20.);
163 
164  fHistPtvsRapidityPromptD0 = new TH2D("fHistPtvsRapidityPromptD0","Prompt D0", 100, 0., 20.,20,-5.,5.);
166  fHistPtvsRapidityPromptLc = new TH2D("fHistPtvsRapidityPromptLc","Prompt Lc", 100, 0., 20.,20,-5.,5.);
168  fHistPtvsRapidityPromptXic0 = new TH2D("fHistPtvsRapidityPromptXic0","Prompt Xic0", 100, 0., 20.,20,-5.,5.);
170  fHistPtvsRapidityFeeddownD0 = new TH2D("fHistPtvsRapidityFeeddownD0","Feeddown D0", 100, 0., 20.,20,-5.,5.);
172  fHistPtvsRapidityFeeddownLc = new TH2D("fHistPtvsRapidityFeeddownLc","Feeddown Lc", 100, 0., 20.,20,-5.,5.);
174  fHistPtvsRapidityFeeddownXic0 = new TH2D("fHistPtvsRapidityFeeddownXic0","Feeddown Xic0", 100, 0., 20.,20,-5.,5.);
176  fHistPtvsRapidityInclusiveD0 = new TH2D("fHistPtvsRapidityInclusiveD0","Inclusive D0", 100, 0., 20.,20,-5.,5.);
178  fHistPtvsRapidityInclusiveLc = new TH2D("fHistPtvsRapidityInclusiveLc","Inclusive Lc", 100, 0., 20.,20,-5.,5.);
180  fHistPtvsRapidityInclusiveXic0 = new TH2D("fHistPtvsRapidityInclusiveXic0","Inclusive Xic0", 100, 0., 20.,20,-5.,5.);
182 
183 }
184 
185 //________________________________________________________________________
187 {
188  //
189  // Main loop
190  // Called for each event
191  //
192 
193  // also a AliEvent...
194  // AliVEvent* mcEvent = MCEvent();
195  AliMCEvent* mcEvent = MCEvent();
196  // AliGenEventHeader * htmp = mcEvent->GenEventHeader();
197 
198  if (!mcEvent) {
199  Printf("ERROR: Could not retrieve MC event");
200  return;
201  }
202 
203  fHistEvt->Fill(0);
204  fHistMult->Fill(mcEvent->GetNumberOfTracks());
205 
206  // Track loop
207  for (Int_t iTrack = 0; iTrack < mcEvent->GetNumberOfTracks(); iTrack++) {
208  AliMCParticle *track = (AliMCParticle*)mcEvent->GetTrack(iTrack);
209  if (!track) {
210  Printf("ERROR: Could not receive track %d", iTrack);
211  continue;
212  }
213  //Bool_t isPrimary = mcEvent->Stack()->IsPhysicalPrimary(iTrack);
214  Int_t pdgcode = abs(track->PdgCode());
215 
216  if(pdgcode==421 || pdgcode==4122 || pdgcode==4132){
217  Int_t pdgarray_chad[100];
218  Int_t labelarray_chad[100];
219  Int_t ngen_chad=0;
220  GetHistory(track,mcEvent,pdgarray_chad,labelarray_chad,ngen_chad);
221  //cout<<isPrimary<<" "<<pdgcode<<" "<<pdgarray_chad[0]<<" "<<pdgarray_chad[1]<<endl;
222 
223  Bool_t isbfd = FromBottom(pdgarray_chad);
224 
225  if(pdgcode==421) fHistPtvsRapidityInclusiveD0->Fill(track->Pt(),track->Y());
226  if(pdgcode==4122) fHistPtvsRapidityInclusiveLc->Fill(track->Pt(),track->Y());
227  if(pdgcode==4132) fHistPtvsRapidityInclusiveXic0->Fill(track->Pt(),track->Y());
228  if(isbfd){
229  if(pdgcode==421) fHistPtvsRapidityFeeddownD0->Fill(track->Pt(),track->Y());
230  if(pdgcode==4122) fHistPtvsRapidityFeeddownLc->Fill(track->Pt(),track->Y());
231  if(pdgcode==4132) fHistPtvsRapidityFeeddownXic0->Fill(track->Pt(),track->Y());
232  }else{
233  if(pdgcode==421) fHistPtvsRapidityPromptD0->Fill(track->Pt(),track->Y());
234  if(pdgcode==4122) fHistPtvsRapidityPromptLc->Fill(track->Pt(),track->Y());
235  if(pdgcode==4132) fHistPtvsRapidityPromptXic0->Fill(track->Pt(),track->Y());
236  }
237 
238  if(TMath::Abs(track->Y())<0.5){
239  if(pdgcode==421) fHistPtInclusiveD0->Fill(track->Pt());
240  if(pdgcode==4122) fHistPtInclusiveLc->Fill(track->Pt());
241  if(pdgcode==4132) fHistPtInclusiveXic0->Fill(track->Pt());
242  if(isbfd){
243  if(pdgcode==421) fHistPtFeeddownD0->Fill(track->Pt());
244  if(pdgcode==4122) fHistPtFeeddownLc->Fill(track->Pt());
245  if(pdgcode==4132) fHistPtFeeddownXic0->Fill(track->Pt());
246  }else{
247  if(pdgcode==421) fHistPtPromptD0->Fill(track->Pt());
248  if(pdgcode==4122) fHistPtPromptLc->Fill(track->Pt());
249  if(pdgcode==4132) fHistPtPromptXic0->Fill(track->Pt());
250  }
251  }
252  }
253  } //track loop
254 
255 
256  // Post output data.
257  PostData(1, fMyOut);
258 }
259 
260 //________________________________________________________________________
262 {
263  //
264  // Draw result to the screen
265  // Called once at the end of the query
266  //
267 
268  fMyOut = dynamic_cast<TList*> (GetOutputData(1));
269 
270  Finalize();
271 }
272 
273 
274 
276 {
277 
278  //
279  // Finalize
280  //
281 
282 }
283 
284 void AliAnalysisTaskCharmBaryonsMC::GetHistory(AliMCParticle *part, AliMCEvent *mcevt, Int_t *pdgarray, Int_t *labelarray, Int_t &ngen)
285 {
286  //
287  // Get history of a particle
288  //
289 
290  for(Int_t i=0;i<100;i++){
291  pdgarray[i] = -9999;
292  labelarray[i] = -9999;
293  }
294  ngen = 0;
295 
296  AliMCParticle *mcprim = part;
297  while(mcprim->GetMother()>=0) {
298  Int_t lab_prim=mcprim->GetMother();
299 
300  AliMCParticle *tmcprim = (AliMCParticle*)mcevt->GetTrack(lab_prim);
301  if(!tmcprim) {
302  break;
303  }
304  if((TMath::Abs(tmcprim->PdgCode())<10) || (TMath::Abs(tmcprim->PdgCode())==21)) break;
305 
306  mcprim = tmcprim;
307 
308  pdgarray[ngen] = mcprim->PdgCode();
309  labelarray[ngen] = lab_prim;
310 
311  ngen ++;
312  if(ngen == 100) break;
313  }
314 }
315 
317  //
318  // Match float
319  //
320  if(TMath::Abs(a-b)<0.00001) return kTRUE;
321  if(TMath::Abs(a+b)<0.00001) return kTRUE;
322  else return kFALSE;
323 }
325  //
326  // Match Int_t
327  //
328  if((a-b) == 0) return kTRUE;
329  if((a+b) == 0) return kTRUE;
330  else return kFALSE;
331 }
332 
334  //
335  // See if the particle is from bottom hadrons
336  //
337  for(int i=0;i<10;i++){
338  if(Match(hitory[i],511)) return kTRUE;
339  if(Match(hitory[i],521)) return kTRUE;
340  if(Match(hitory[i],531)) return kTRUE;
341  if(Match(hitory[i],5122)) return kTRUE;
342  if(Match(hitory[i],5232)) return kTRUE;
343  if(Match(hitory[i],5132)) return kTRUE;
344  if(Match(hitory[i],5332)) return kTRUE;
345  }
346  return kFALSE;
347 }
TH2D * fHistPtvsRapidityFeeddownLc
! Lc pT distribution in |y|<0.5
TH1D * fHistPtFeeddownXic0
! Xic0 pT distribution in |y|<0.5
ClassImp(AliAnalysisTaskCharmBaryonsMC) AliAnalysisTaskCharmBaryonsMC
void GetHistory(AliMCParticle *part, AliMCEvent *mcevt, Int_t *pdgarray, Int_t *labelarray, Int_t &ngen)
TH1D * fHistPtFeeddownLc
! Lc pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityInclusiveXic0
! Xic0 pT distribution in |y|<0.5
int Int_t
Definition: External.C:63
TH2D * fHistPtvsRapidityPromptXic0
! Xic0 pT distribution in |y|<0.5
TH1D * fHistPtPromptLc
! Lc pT distribution in |y|<0.5
float Float_t
Definition: External.C:68
Definition: External.C:228
Definition: External.C:212
virtual void UserExec(Option_t *option)
TH1D * fHistPtFeeddownD0
! D0 pT distribution in |y|<0.5
TH1D * fHistPtPromptD0
! D0 pT distribution in |y|<0.5
TH1D * fHistPtPromptXic0
! Xic0 pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityPromptLc
! Lc pT distribution in |y|<0.5
TH1D * fHistPtInclusiveXic0
! Xic0 pT distribution in |y|<0.5
TH1D * fHistPtInclusiveLc
! Lc pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityFeeddownD0
! D0 pT distribution in |y|<0.5
const char Option_t
Definition: External.C:48
TH1D * fHistPtInclusiveD0
! D0 pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityInclusiveLc
! Lc pT distribution in |y|<0.5
bool Bool_t
Definition: External.C:53
TH2D * fHistPtvsRapidityFeeddownXic0
! Xic0 pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityInclusiveD0
! D0 pT distribution in |y|<0.5
TH2D * fHistPtvsRapidityPromptD0
! D0 pT distribution in |y|<0.5
TList * fMyOut
! list of output histos