1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //_________________________________________________________________________
17 // Do photon/pi0 analysis for isolation and correlation
18 // at the generator level. Only for kine stack (ESDs)
21 // -- Author: Gustavo Conesa (LPSC-CNRS-Grenoble)
22 //////////////////////////////////////////////////////////////////////////////
24 // --- ROOT system ---
26 #include "TParticle.h"
27 #include "TDatabasePDG.h"
29 //---- ANALYSIS system ----
30 #include "AliAnaGeneratorKine.h"
32 #include "AliGenPythiaEventHeader.h"
34 ClassImp(AliAnaGeneratorKine)
37 //__________________________________________
38 AliAnaGeneratorKine::AliAnaGeneratorKine() :
39 AliAnaCaloTrackCorrBaseClass(),
41 fParton2(0), fParton3(0),
42 fParton6(0), fParton7(0),
45 fhPtHard(0), fhPtParton(0), fhPtJet(0),
46 fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0),
47 fhPtPhoton(0), fhPtPi0(0)
51 //Initialize parameters
54 for(Int_t i = 0; i < 4; i++)
56 fhPtPhotonLeading[i] = fhPtPi0Leading[i] = 0;
57 fhPtPhotonLeadingIsolated[i] = fhPtPi0LeadingIsolated[i] = 0;
58 fhZHardPhotonLeading[i] = fhZHardPi0Leading[i] = 0;
59 fhZHardPhotonLeadingIsolated[i] = fhZHardPi0LeadingIsolated[i] = 0;
60 fhZPartonPhotonLeading[i] = fhZPartonPi0Leading[i] = 0;
61 fhZPartonPhotonLeadingIsolated[i] = fhZPartonPi0LeadingIsolated[i] = 0;
62 fhZJetPhotonLeading[i] = fhZJetPi0Leading[i] = 0;
63 fhZJetPhotonLeadingIsolated[i] = fhZJetPi0LeadingIsolated[i] = 0;
64 fhXEPhotonLeading[i] = fhXEPi0Leading[i] = 0;
65 fhXEPhotonLeadingIsolated[i] = fhXEPi0LeadingIsolated[i] = 0;
66 fhXEUEPhotonLeading[i] = fhXEUEPi0Leading[i] = 0;
67 fhXEUEPhotonLeadingIsolated[i] = fhXEUEPi0LeadingIsolated[i] = 0;
69 fhPtPartonTypeNearPhotonLeading[i] = fhPtPartonTypeNearPi0Leading[i] = 0;
70 fhPtPartonTypeNearPhotonLeadingIsolated[i] = fhPtPartonTypeNearPi0LeadingIsolated[i] = 0;
72 fhPtPartonTypeAwayPhotonLeading[i] = fhPtPartonTypeAwayPi0Leading[i] = 0;
73 fhPtPartonTypeAwayPhotonLeadingIsolated[i] = fhPtPartonTypeAwayPi0LeadingIsolated[i] = 0;
79 //_______________________________________________________________________________
80 Bool_t AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger,
81 const Int_t indexTrig,
83 const Bool_t leading[4],
84 const Bool_t isolated[4],
87 //Correlate trigger with partons or jets, get z
89 //Get the index of the mother
90 iparton = (fStack->Particle(indexTrig))->GetFirstMother();
91 TParticle * mother = fStack->Particle(iparton);
94 iparton = mother->GetFirstMother();
95 if(iparton < 0) { printf("AliAnaGeneratorKine::CorrelateWithPartonOrJet() - Negative index, skip event\n"); return kFALSE; }
96 mother = fStack->Particle(iparton);
99 //printf("Mother is parton %d with pdg %d\n",imom,mother->GetPdgCode());
103 //printf("This particle is not from hard process - pdg %d, parton index %d\n",pdgTrig, iparton);
107 Float_t ptTrig = trigger.Pt();
108 Float_t partonPt = fParton6->Pt();
109 Float_t jetPt = fJet6.Pt();
112 partonPt = fParton6->Pt();
116 //Get id of parton in near and away side
123 //printf("parton 6 pdg = %d, parton 7 pdg = %d\n",fParton6->GetPdgCode(),fParton7->GetPdgCode());
127 nearPDG = fParton6->GetPdgCode();
128 awayPDG = fParton7->GetPdgCode();
132 nearPDG = fParton7->GetPdgCode();
133 awayPDG = fParton6->GetPdgCode();
136 if (nearPDG == 22) near = 0;
137 else if(nearPDG == 21) near = 1;
140 if (awayPDG == 22) away = 0;
141 else if(awayPDG == 21) away = 1;
144 for( Int_t i = 0; i < 4; i++ )
150 fhPtPartonTypeNearPi0Leading[i]->Fill(ptTrig,near);
151 fhPtPartonTypeAwayPi0Leading[i]->Fill(ptTrig,away);
154 fhPtPartonTypeNearPi0LeadingIsolated[i]->Fill(ptTrig,near);
155 fhPtPartonTypeAwayPi0LeadingIsolated[i]->Fill(ptTrig,away);
163 fhPtPartonTypeNearPhotonLeading[i]->Fill(ptTrig,near);
164 fhPtPartonTypeAwayPhotonLeading[i]->Fill(ptTrig,away);
167 fhPtPartonTypeNearPhotonLeadingIsolated[i]->Fill(ptTrig,near);
168 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->Fill(ptTrig,away);
177 fhPtPartonPtHard->Fill(fPtHard, partonPt/fPtHard);
178 fhPtJetPtHard ->Fill(fPtHard, jetPt/fPtHard);
179 fhPtJetPtParton ->Fill(fPtHard, jetPt/partonPt);
181 Float_t zHard = ptTrig / fPtHard;
182 Float_t zPart = ptTrig / partonPt;
183 Float_t zJet = ptTrig / jetPt;
185 //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard);
187 //printf("Z : hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
189 for( Int_t i = 0; i < 4; i++ )
195 fhZHardPi0Leading[i] ->Fill(ptTrig,zHard);
196 fhZPartonPi0Leading[i]->Fill(ptTrig,zPart);
197 fhZJetPi0Leading[i] ->Fill(ptTrig,zJet );
201 fhZHardPi0LeadingIsolated[i] ->Fill(ptTrig,zHard);
202 fhZPartonPi0LeadingIsolated[i]->Fill(ptTrig,zPart);
203 fhZJetPi0LeadingIsolated[i] ->Fill(ptTrig,zJet);
211 fhZHardPhotonLeading[i] ->Fill(ptTrig,zHard);
212 fhZPartonPhotonLeading[i]->Fill(ptTrig,zPart);
213 fhZJetPhotonLeading[i] ->Fill(ptTrig,zJet );
217 fhZHardPhotonLeadingIsolated[i] ->Fill(ptTrig,zHard);
218 fhZPartonPhotonLeadingIsolated[i]->Fill(ptTrig,zPart);
219 fhZJetPhotonLeadingIsolated[i] ->Fill(ptTrig,zJet);
229 //____________________________________________________
230 TList * AliAnaGeneratorKine::GetCreateOutputObjects()
232 // Create histograms to be saved in output file
234 TList * outputContainer = new TList() ;
235 outputContainer->SetName("GenKineHistos") ;
237 Int_t nptbins = GetHistogramRanges()->GetHistoPtBins();
238 Float_t ptmax = GetHistogramRanges()->GetHistoPtMax();
239 Float_t ptmin = GetHistogramRanges()->GetHistoPtMin();
241 fhPtHard = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax);
242 fhPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
243 outputContainer->Add(fhPtHard);
245 fhPtParton = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax);
246 fhPtParton->SetXTitle("p_{T}^{parton} (GeV/c)");
247 outputContainer->Add(fhPtParton);
249 fhPtJet = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax);
250 fhPtJet->SetXTitle("p_{T}^{jet} (GeV/c)");
251 outputContainer->Add(fhPtJet);
253 fhPtPartonPtHard = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
254 fhPtPartonPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
255 fhPtPartonPtHard->SetYTitle("p_{T}^{parton}/p_{T}^{hard}");
256 outputContainer->Add(fhPtPartonPtHard);
258 fhPtJetPtHard = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
259 fhPtJetPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
260 fhPtJetPtHard->SetYTitle("p_{T}^{jet}/p_{T}^{hard}");
261 outputContainer->Add(fhPtJetPtHard);
263 fhPtJetPtParton = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2);
264 fhPtJetPtParton->SetXTitle("p_{T}^{hard} (GeV/c)");
265 fhPtJetPtParton->SetYTitle("p_{T}^{jet}/p_{T}^{parton}");
266 outputContainer->Add(fhPtJetPtParton);
269 fhPtPhoton = new TH1F("hPtPhoton","Input Photon",nptbins,ptmin,ptmax);
270 fhPtPhoton->SetXTitle("p_{T} (GeV/c)");
271 outputContainer->Add(fhPtPhoton);
273 fhPtPi0 = new TH1F("hPtPi0","Input Pi0",nptbins,ptmin,ptmax);
274 fhPtPi0->SetXTitle("p_{T} (GeV/c)");
275 outputContainer->Add(fhPtPi0);
277 TString name[] = {"","_EMC","_Photon","_EMC_Photon"};
278 TString title[] = {"",", neutral in EMCal",", neutral only photon like",", neutral in EMCal and only photon like"};
280 for(Int_t i = 0; i < 4; i++)
285 fhPtPhotonLeading[i] = new TH1F(Form("hPtPhotonLeading%s",name[i].Data()),
286 Form("Photon : Leading of all particles%s",title[i].Data()),
287 nptbins,ptmin,ptmax);
288 fhPtPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
289 outputContainer->Add(fhPtPhotonLeading[i]);
291 fhPtPi0Leading[i] = new TH1F(Form("hPtPi0Leading%s",name[i].Data()),
292 Form("Pi0 : Leading of all particles%s",title[i].Data()),
293 nptbins,ptmin,ptmax);
294 fhPtPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
295 outputContainer->Add(fhPtPi0Leading[i]);
297 fhPtPhotonLeadingIsolated[i] = new TH1F(Form("hPtPhotonLeadingIsolated%s",name[i].Data()),
298 Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
299 nptbins,ptmin,ptmax);
300 fhPtPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
301 outputContainer->Add(fhPtPhotonLeadingIsolated[i]);
303 fhPtPi0LeadingIsolated[i] = new TH1F(Form("hPtPi0LeadingIsolated%s",name[i].Data()),
304 Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
305 nptbins,ptmin,ptmax);
306 fhPtPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
307 outputContainer->Add(fhPtPi0LeadingIsolated[i]);
311 fhPtPartonTypeNearPhotonLeading[i] = new TH2F(Form("hPtPartonTypeNearPhotonLeading%s",name[i].Data()),
312 Form("Photon : Leading of all particles%s",title[i].Data()),
313 nptbins,ptmin,ptmax,3,0,3);
314 fhPtPartonTypeNearPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
315 fhPtPartonTypeNearPhotonLeading[i]->SetYTitle("Parton type");
316 fhPtPartonTypeNearPhotonLeading[i]->GetYaxis()->SetBinLabel(1,"#gamma");
317 fhPtPartonTypeNearPhotonLeading[i]->GetYaxis()->SetBinLabel(2,"g");
318 fhPtPartonTypeNearPhotonLeading[i]->GetYaxis()->SetBinLabel(3,"q");
319 outputContainer->Add(fhPtPartonTypeNearPhotonLeading[i]);
321 fhPtPartonTypeNearPi0Leading[i] = new TH2F(Form("hPtPartonTypeNearPi0Leading%s",name[i].Data()),
322 Form("Pi0 : Leading of all particles%s",title[i].Data()),
323 nptbins,ptmin,ptmax,3,0,3);
324 fhPtPartonTypeNearPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
325 fhPtPartonTypeNearPi0Leading[i]->SetYTitle("Parton type");
326 fhPtPartonTypeNearPi0Leading[i]->GetYaxis()->SetBinLabel(1,"#gamma");
327 fhPtPartonTypeNearPi0Leading[i]->GetYaxis()->SetBinLabel(2,"g");
328 fhPtPartonTypeNearPi0Leading[i]->GetYaxis()->SetBinLabel(3,"q");
329 outputContainer->Add(fhPtPartonTypeNearPi0Leading[i]);
331 fhPtPartonTypeNearPhotonLeadingIsolated[i] = new TH2F(Form("hPtPartonTypeNearPhotonLeadingIsolated%s",name[i].Data()),
332 Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
333 nptbins,ptmin,ptmax,3,0,3);
334 fhPtPartonTypeNearPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
335 fhPtPartonTypeNearPhotonLeadingIsolated[i]->SetYTitle("Parton type");
336 fhPtPartonTypeNearPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(1,"#gamma");
337 fhPtPartonTypeNearPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(2,"g");
338 fhPtPartonTypeNearPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(3,"q");
339 outputContainer->Add(fhPtPartonTypeNearPhotonLeadingIsolated[i]);
341 fhPtPartonTypeNearPi0LeadingIsolated[i] = new TH2F(Form("hPtPartonTypeNearPi0LeadingIsolated%s",name[i].Data()),
342 Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
343 nptbins,ptmin,ptmax,3,0,3);
344 fhPtPartonTypeNearPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
345 fhPtPartonTypeNearPi0LeadingIsolated[i]->SetYTitle("Parton type");
346 fhPtPartonTypeNearPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(1,"#gamma");
347 fhPtPartonTypeNearPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(2,"g");
348 fhPtPartonTypeNearPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(3,"q");
349 outputContainer->Add(fhPtPartonTypeNearPi0LeadingIsolated[i]);
354 fhPtPartonTypeAwayPhotonLeading[i] = new TH2F(Form("hPtPartonTypeAwayPhotonLeading%s",name[i].Data()),
355 Form("Photon : Leading of all particles%s",title[i].Data()),
356 nptbins,ptmin,ptmax,3,0,3);
357 fhPtPartonTypeAwayPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
358 fhPtPartonTypeAwayPhotonLeading[i]->SetYTitle("Parton type");
359 fhPtPartonTypeAwayPhotonLeading[i]->GetYaxis()->SetBinLabel(1,"#gamma");
360 fhPtPartonTypeAwayPhotonLeading[i]->GetYaxis()->SetBinLabel(2,"g");
361 fhPtPartonTypeAwayPhotonLeading[i]->GetYaxis()->SetBinLabel(3,"q");
362 outputContainer->Add(fhPtPartonTypeAwayPhotonLeading[i]);
364 fhPtPartonTypeAwayPi0Leading[i] = new TH2F(Form("hPtPartonTypeAwayPi0Leading%s",name[i].Data()),
365 Form("Pi0 : Leading of all particles%s",title[i].Data()),
366 nptbins,ptmin,ptmax,3,0,3);
367 fhPtPartonTypeAwayPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
368 fhPtPartonTypeAwayPi0Leading[i]->SetYTitle("Parton type");
369 fhPtPartonTypeAwayPi0Leading[i]->GetYaxis()->SetBinLabel(1,"#gamma");
370 fhPtPartonTypeAwayPi0Leading[i]->GetYaxis()->SetBinLabel(2,"g");
371 fhPtPartonTypeAwayPi0Leading[i]->GetYaxis()->SetBinLabel(3,"q");
372 outputContainer->Add(fhPtPartonTypeAwayPi0Leading[i]);
374 fhPtPartonTypeAwayPhotonLeadingIsolated[i] = new TH2F(Form("hPtPartonTypeAwayPhotonLeadingIsolated%s",name[i].Data()),
375 Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
376 nptbins,ptmin,ptmax,3,0,3);
377 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
378 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->SetYTitle("Parton type");
379 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(1,"#gamma");
380 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(2,"g");
381 fhPtPartonTypeAwayPhotonLeadingIsolated[i]->GetYaxis()->SetBinLabel(3,"q");
382 outputContainer->Add(fhPtPartonTypeAwayPhotonLeadingIsolated[i]);
384 fhPtPartonTypeAwayPi0LeadingIsolated[i] = new TH2F(Form("hPtPartonTypeAwayPi0LeadingIsolated%s",name[i].Data()),
385 Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
386 nptbins,ptmin,ptmax,3,0,3);
387 fhPtPartonTypeAwayPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
388 fhPtPartonTypeAwayPi0LeadingIsolated[i]->SetYTitle("Parton type");
389 fhPtPartonTypeAwayPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(1,"#gamma");
390 fhPtPartonTypeAwayPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(2,"g");
391 fhPtPartonTypeAwayPi0LeadingIsolated[i]->GetYaxis()->SetBinLabel(3,"q");
392 outputContainer->Add(fhPtPartonTypeAwayPi0LeadingIsolated[i]);
396 fhZHardPhotonLeading[i] = new TH2F(Form("hZHardPhotonLeading%s",name[i].Data()),
397 Form("Z-Hard of Photon : Leading of all particles%s",title[i].Data()),
398 nptbins,ptmin,ptmax,200,0,2);
399 fhZHardPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
400 fhZHardPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
401 outputContainer->Add(fhZHardPhotonLeading[i]);
403 fhZHardPi0Leading[i] = new TH2F(Form("hZHardPi0Leading%s",name[i].Data()),
404 Form("Z-Hard of Pi0 : Leading of all particles%s",title[i].Data()),
405 nptbins,ptmin,ptmax,200,0,2);
406 fhZHardPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
407 fhZHardPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
408 outputContainer->Add(fhZHardPi0Leading[i]);
410 fhZHardPhotonLeadingIsolated[i] = new TH2F(Form("hZHardPhotonLeadingIsolated%s",name[i].Data()),
411 Form("Z-Hard of Photon : Leading of all particles%s, isolated",title[i].Data()),
412 nptbins,ptmin,ptmax,200,0,2);
413 fhZHardPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
414 fhZHardPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
415 outputContainer->Add(fhZHardPhotonLeadingIsolated[i]);
417 fhZHardPi0LeadingIsolated[i] = new TH2F(Form("hZHardPi0LeadingIsolated%s",name[i].Data()),
418 Form("Z-Hard of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
419 nptbins,ptmin,ptmax,200,0,2);
420 fhZHardPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
421 fhZHardPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
422 outputContainer->Add(fhZHardPi0LeadingIsolated[i]);
426 fhZPartonPhotonLeading[i] = new TH2F(Form("hZPartonPhotonLeading%s",name[i].Data()),
427 Form("Z-Parton of Photon : Leading of all particles%s",title[i].Data()),
428 nptbins,ptmin,ptmax,200,0,2);
429 fhZPartonPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
430 fhZPartonPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
431 outputContainer->Add(fhZPartonPhotonLeading[i]);
433 fhZPartonPi0Leading[i] = new TH2F(Form("hZPartonPi0Leading%s",name[i].Data()),
434 Form("Z-Parton of Pi0 : Leading of all particles%s",title[i].Data()),
435 nptbins,ptmin,ptmax,200,0,2);
436 fhZPartonPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
437 fhZPartonPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
438 outputContainer->Add(fhZPartonPi0Leading[i]);
440 fhZPartonPhotonLeadingIsolated[i] = new TH2F(Form("hZPartonPhotonLeadingIsolated%s",name[i].Data()),
441 Form("Z-Parton of Photon : Leading of all particles%s, isolated",title[i].Data()),
442 nptbins,ptmin,ptmax,200,0,2);
443 fhZPartonPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
444 fhZPartonPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
445 outputContainer->Add(fhZPartonPhotonLeadingIsolated[i]);
447 fhZPartonPi0LeadingIsolated[i] = new TH2F(Form("hZPartonPi0LeadingIsolated%s",name[i].Data()),
448 Form("Z-Parton of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
449 nptbins,ptmin,ptmax,200,0,2);
450 fhZPartonPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
451 fhZPartonPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
452 outputContainer->Add(fhZPartonPi0LeadingIsolated[i]);
457 fhZJetPhotonLeading[i] = new TH2F(Form("hZJetPhotonLeading%s",name[i].Data()),
458 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
459 nptbins,ptmin,ptmax,200,0,2);
460 fhZJetPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
461 fhZJetPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
462 outputContainer->Add(fhZJetPhotonLeading[i]);
464 fhZJetPi0Leading[i] = new TH2F(Form("hZJetPi0Leading%s",name[i].Data()),
465 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
466 nptbins,ptmin,ptmax,200,0,2);
467 fhZJetPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
468 fhZJetPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
469 outputContainer->Add(fhZJetPi0Leading[i]);
471 fhZJetPhotonLeadingIsolated[i] = new TH2F(Form("hZJetPhotonLeadingIsolated%s",name[i].Data()),
472 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
473 nptbins,ptmin,ptmax,200,0,2);
474 fhZJetPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
475 fhZJetPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
476 outputContainer->Add(fhZJetPhotonLeadingIsolated[i]);
478 fhZJetPi0LeadingIsolated[i] = new TH2F(Form("hZJetPi0LeadingIsolated%s",name[i].Data()),
479 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
480 nptbins,ptmin,ptmax,200,0,2);
481 fhZJetPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
482 fhZJetPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
483 outputContainer->Add(fhZJetPi0LeadingIsolated[i]);
488 fhXEPhotonLeading[i] = new TH2F(Form("hXEPhotonLeading%s",name[i].Data()),
489 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
490 nptbins,ptmin,ptmax,200,0,2);
491 fhXEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
492 fhXEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
493 outputContainer->Add(fhXEPhotonLeading[i]);
495 fhXEPi0Leading[i] = new TH2F(Form("hXEPi0Leading%s",name[i].Data()),
496 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
497 nptbins,ptmin,ptmax,200,0,2);
498 fhXEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
499 fhXEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
500 outputContainer->Add(fhXEPi0Leading[i]);
502 fhXEPhotonLeadingIsolated[i] = new TH2F(Form("hXEPhotonLeadingIsolated%s",name[i].Data()),
503 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
504 nptbins,ptmin,ptmax,200,0,2);
505 fhXEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
506 fhXEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
507 outputContainer->Add(fhXEPhotonLeadingIsolated[i]);
509 fhXEPi0LeadingIsolated[i] = new TH2F(Form("hXEPi0LeadingIsolated%s",name[i].Data()),
510 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
511 nptbins,ptmin,ptmax,200,0,2);
512 fhXEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
513 fhXEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
514 outputContainer->Add(fhXEPi0LeadingIsolated[i]);
519 fhXEUEPhotonLeading[i] = new TH2F(Form("hXEUEPhotonLeading%s",name[i].Data()),
520 Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
521 nptbins,ptmin,ptmax,200,0,2);
522 fhXEUEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
523 fhXEUEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
524 outputContainer->Add(fhXEUEPhotonLeading[i]);
526 fhXEUEPi0Leading[i] = new TH2F(Form("hXEUEPi0Leading%s",name[i].Data()),
527 Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
528 nptbins,ptmin,ptmax,200,0,2);
529 fhXEUEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
530 fhXEUEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
531 outputContainer->Add(fhXEUEPi0Leading[i]);
533 fhXEUEPhotonLeadingIsolated[i] = new TH2F(Form("hXEUEPhotonLeadingIsolated%s",name[i].Data()),
534 Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
535 nptbins,ptmin,ptmax,200,0,2);
536 fhXEUEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
537 fhXEUEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
538 outputContainer->Add(fhXEUEPhotonLeadingIsolated[i]);
540 fhXEUEPi0LeadingIsolated[i] = new TH2F(Form("hXEUEPi0LeadingIsolated%s",name[i].Data()),
541 Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
542 nptbins,ptmin,ptmax,200,0,2);
543 fhXEUEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
544 fhXEUEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
545 outputContainer->Add(fhXEUEPi0LeadingIsolated[i]);
549 return outputContainer;
553 //____________________________________________
554 void AliAnaGeneratorKine::GetPartonsAndJets()
556 // Fill data members with partons,jets and generated pt hard
558 fStack = GetMCStack() ;
562 printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - No Stack available, STOP\n");
566 fParton2 = fStack->Particle(2) ;
567 fParton3 = fStack->Particle(3) ;
568 fParton6 = fStack->Particle(6) ;
569 fParton7 = fStack->Particle(7) ;
571 Float_t p6phi = fParton6->Phi();
572 if(p6phi < 0) p6phi +=TMath::TwoPi();
573 Float_t p7phi = fParton7->Phi();
574 if(p7phi < 0) p7phi +=TMath::TwoPi();
576 //printf("parton6: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton6->Pt(),fParton6->Eta(),p6phi, fParton6->GetPdgCode());
577 //printf("parton7: pt %2.2f, eta %2.2f, phi %2.2f with pdg %d\n",fParton7->Pt(),fParton7->Eta(),p7phi, fParton7->GetPdgCode());
579 // Get the jets, only for pythia
580 if(!strcmp(GetReader()->GetGenEventHeader()->ClassName(), "AliGenPythiaEventHeader"))
582 AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) GetReader()->GetGenEventHeader();
584 fPtHard = pygeh->GetPtHard();
586 //printf("pt Hard %2.2f\n",fPtHard);
588 const Int_t nTriggerJets = pygeh->NTriggerJets();
590 Float_t tmpjet[]={0,0,0,0};
592 // select the closest jet to parton
596 for(Int_t ijet = 0; ijet< nTriggerJets; ijet++)
598 pygeh->TriggerJet(ijet, tmpjet);
600 TLorentzVector jet(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
601 Float_t jphi = jet.Phi();
602 if(jphi < 0) jphi +=TMath::TwoPi();
604 Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, jet.Eta() , jphi) ;
605 Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, jet.Eta() , jphi) ;
607 //printf("jet %d: pt %2.2f, eta %2.2f, phi %2.2f, r6 %2.2f, r7 %2.2f\n",ijet,jet.Pt(),jet.Eta(),jphi,radius6, radius7);
623 //printf("jet6: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet6.Pt(),fJet6.Eta(),fJet6.Phi());
624 //printf("jet7: pt %2.2f, eta %2.2f, phi %2.2f\n",fJet7.Pt(),fJet7.Eta(),fJet6.Phi());
628 fhPtHard ->Fill(fPtHard);
629 fhPtJet ->Fill(fJet6.Pt());
630 fhPtJet ->Fill(fJet7.Pt());
631 fhPtParton ->Fill(fParton6->Pt());
632 fhPtParton ->Fill(fParton7->Pt());
636 //___________________________________________________________
637 void AliAnaGeneratorKine::GetXE(const TLorentzVector trigger,
638 const Int_t indexTrig,
640 const Bool_t leading[4],
641 const Bool_t isolated[4],
645 // Calculate the real XE and the UE XE
647 Float_t ptThresTrack = 0.2;
649 Float_t ptTrig = trigger.Pt();
650 Float_t etaTrig = trigger.Eta();
651 Float_t phiTrig = trigger.Phi();
652 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
654 //Loop on primaries, start from position 8, no partons
655 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
657 TParticle * particle = fStack->Particle(ipr) ;
659 if(ipr==indexTrig) continue;
662 Int_t pdg = particle->GetPdgCode();
663 Int_t status = particle->GetStatusCode();
665 // Compare trigger with final state particles
666 if( status != 1) continue ;
668 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
670 if(charge==0) continue; // construct xe only with charged
672 Float_t pt = particle->Pt();
673 Float_t eta = particle->Eta();
674 Float_t phi = particle->Phi();
675 if(phi < 0 ) phi += TMath::TwoPi();
677 if( pt < ptThresTrack) continue ;
679 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
682 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
684 Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
686 //Get the index of the mother
687 Int_t ipartonAway = particle->GetFirstMother();
688 if(ipartonAway < 0) return;
689 TParticle * mother = fStack->Particle(ipartonAway);
690 while (ipartonAway > 7)
692 ipartonAway = mother->GetFirstMother();
693 if(ipartonAway < 0) break;
694 mother = fStack->Particle(ipartonAway);
697 if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway)
699 //printf("xE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
700 if(radius > 1 ) continue; // avoid particles too far from trigger
702 for( Int_t i = 0; i < 4; i++ )
708 fhXEPi0Leading[i] ->Fill(ptTrig,xe);
712 fhXEPi0LeadingIsolated[i] ->Fill(ptTrig,xe);
720 fhXEPhotonLeading[i] ->Fill(ptTrig,xe);
724 fhXEPhotonLeadingIsolated[i] ->Fill(ptTrig,xe);
731 if(ipartonAway!=6 && ipartonAway!=7)
734 //printf("xE UE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
736 for( Int_t i = 0; i < 4; i++ )
742 fhXEUEPi0Leading[i] ->Fill(ptTrig,xe);
746 fhXEUEPi0LeadingIsolated[i] ->Fill(ptTrig,xe);
754 fhXEUEPhotonLeading[i] ->Fill(ptTrig,xe);
758 fhXEUEPhotonLeadingIsolated[i] ->Fill(ptTrig,xe);
770 //________________________________________
771 void AliAnaGeneratorKine::InitParameters()
774 //Initialize the parameters of the analysis.
775 AddToHistogramsName("AnaGenKine_");
779 //___________________________________________________________________________
780 void AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
781 const Int_t indexTrig,
786 // Check if the trigger is the leading particle
787 // In case of neutral particles check all neutral or neutral in EMCAL acceptance
789 Float_t ptMaxCharged = 0; // all charged
790 Float_t ptMaxNeutral = 0; // all neutral
791 Float_t ptMaxNeutEMCAL = 0; // for neutral, select them in EMCAL acceptance
792 Float_t ptMaxNeutPhot = 0; // for neutral, take only photons
793 Float_t ptMaxNeutEMCALPhot = 0; // for neutral, take only photons in EMCAL acceptance
805 Float_t ptTrig = trigger.Pt();
806 Float_t etaTrig = trigger.Eta();
807 Float_t phiTrig = trigger.Phi();
808 if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
810 // Minimum track or cluster energy
811 Float_t ptThresTrack = 0.2;
812 Float_t ptThresCalo = 0.3;
815 Float_t ptThresIC = 0.5;
816 Float_t rThresIC = 0.4;
818 Int_t nICNeutral = 0;
819 Int_t nICNeutEMCAL = 0;
820 Int_t nICNeutPhot = 0;
821 Int_t nICNeutEMCALPhot = 0;
823 //Loop on primaries, start from position 8, no partons
824 for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
826 if(ipr == indexTrig) continue;
827 TParticle * particle = fStack->Particle(ipr) ;
829 Int_t imother= particle->GetFirstMother();
830 //printf("Leading ipr %d - mother %d\n",ipr, imother);
832 if(imother==indexTrig) continue ;
834 Int_t pdg = particle->GetPdgCode();
835 Int_t status = particle->GetStatusCode();
837 // Compare trigger with final state particles
838 if( status != 1) continue ;
840 Float_t pt = particle->Pt();
841 Float_t eta = particle->Eta();
842 Float_t phi = particle->Phi();
843 if(phi < 0 ) phi += TMath::TwoPi();
845 if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
847 Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
850 Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
854 if(pt < ptThresCalo) continue ;
856 if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
857 if( pt > ptThresIC && radius < rThresIC ) nICNeutral++ ;
859 Bool_t phPDG = kFALSE;
860 if(pdg==22 || pdg==111) phPDG = kTRUE;
862 //if(pt > ptTrig) printf(" --- pdg %d, phPDG %d pT %2.2f, pTtrig %2.2f, eta %2.2f, phi %2.2f\n",pdg,phPDG,pt,ptTrig,particle->Eta(), particle->Phi()*TMath::RadToDeg());
865 if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
866 if( pt > ptThresIC && radius < rThresIC ) nICNeutPhot++ ;
870 Bool_t inEMCAL = kTRUE;
871 if(TMath::Abs(particle->Eta()) > 0.7
872 || particle->Phi() > TMath::DegToRad()*180
873 || particle->Phi() < TMath::DegToRad()*80 ) inEMCAL = kFALSE ;
877 if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
878 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCAL++ ;
882 if( ptMaxNeutEMCALPhot < pt ) ptMaxNeutEMCALPhot = pt;
883 if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCALPhot++ ;
889 if( pt < ptThresTrack) continue ;
891 if( ptMaxCharged < pt ) ptMaxCharged = pt;
893 if( pt > ptThresIC && radius < rThresIC )
895 //printf("UE track? pTtrig %2.2f, pt %2.2f, etaTrig %2.2f, eta %2.2f, phiTrig %2.2f, phi %2.2f, radius %2.2f\n",
896 // ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
904 if(ptTrig > ptMaxCharged)
906 //printf("pt charged %2.2f, pt neutral %2.2f, pt neutral emcal %2.2f, pt photon %2.2f, pt photon emcal %2.2f\n",
907 // ptMaxCharged, ptMaxNeutral, ptMaxNeutEMCAL,ptMaxNeutPhot, ptMaxNeutEMCALPhot);
908 if(ptTrig > ptMaxNeutral ) leading[0] = kTRUE ;
909 if(ptTrig > ptMaxNeutEMCAL ) leading[1] = kTRUE ;
910 if(ptTrig > ptMaxNeutPhot ) leading[2] = kTRUE ;
911 if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
914 //printf("N in cone over threshold : tracks %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n",
915 // nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
917 // Isolation decision
920 if(nICNeutral == 0 ) isolated[0] = kTRUE ;
921 if(nICNeutEMCAL == 0 ) isolated[1] = kTRUE ;
922 if(nICNeutPhot == 0 ) isolated[2] = kTRUE ;
923 if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
926 // Fill histograms if conditions apply for all 4 cases
927 for( Int_t i = 0; i < 4; i++ )
933 fhPtPi0Leading[i]->Fill(ptTrig);
934 if(isolated[i]) fhPtPi0LeadingIsolated[i]->Fill(ptTrig);
941 fhPtPhotonLeading[i]->Fill(ptTrig);
942 if(isolated[i]) fhPtPhotonLeadingIsolated[i]->Fill(ptTrig);
949 //_____________________________________________________
950 void AliAnaGeneratorKine::MakeAnalysisFillHistograms()
952 //Particle-Parton Correlation Analysis, fill histograms
954 TLorentzVector trigger;
958 for(Int_t ipr = 0; ipr < fStack->GetNprimary(); ipr ++ )
960 TParticle * particle = fStack->Particle(ipr) ;
962 Int_t pdgTrig = particle->GetPdgCode();
963 Int_t statusTrig = particle->GetStatusCode();
964 Int_t imother = particle->GetFirstMother();
965 Float_t ptTrig = particle->Pt();
967 // Select final state photons (prompt, fragmentation) or pi0s
969 //Check the origin of the photon, accept if prompt or initial/final state radiation
970 if(pdgTrig == 22 && statusTrig == 1)
972 if(imother > 8) continue;
974 // If not photon, trigger on pi0
975 else if(pdgTrig != 111) continue;
977 // Acceptance and kinematical cuts
978 if( ptTrig < GetMinPt() ) continue ;
980 //EMCAL acceptance, a bit less
981 if(TMath::Abs(particle->Eta()) > 0.6) continue ;
982 if(particle->Phi() > TMath::DegToRad()*176) continue ;
983 if(particle->Phi() < TMath::DegToRad()*74 ) continue ;
985 // printf("Particle %d : pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
986 // ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
990 // printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
994 if (pdgTrig==22 ) fhPtPhoton->Fill(ptTrig);
995 else if(pdgTrig==111) fhPtPi0 ->Fill(ptTrig);
997 // Check if it is leading
1001 particle->Momentum(trigger);
1003 IsLeadingAndIsolated(trigger, ipr, pdgTrig, leading, isolated);
1006 Int_t ok = CorrelateWithPartonOrJet(trigger, ipr, pdgTrig, leading, isolated, iparton);
1009 GetXE(trigger,ipr,pdgTrig,leading,isolated,iparton) ;
1013 if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - End fill histograms \n");