AliPhysics  fde8a9f (fde8a9f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetExtractor.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: R. Haake. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
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  **************************************************************************/
15 
16 #include <cstdlib>
17 #include <algorithm>
18 #include <vector>
19 #include <TClonesArray.h>
20 #include <TF1.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <THn.h>
25 #include <TTree.h>
26 #include <TList.h>
27 #include <TLorentzVector.h>
28 
29 #include "AliEmcalPythiaInfo.h"
30 #include "AliMCEvent.h"
31 #include "AliPythia.h"
32 #include "AliStack.h"
33 
34 #include "TObjArray.h"
35 #include "AliESDVertex.h"
36 #include "AliAODVertex.h"
37 #include "AliAODPid.h"
38 
39 #include "AliVTrack.h"
40 #include "AliVHeader.h"
41 #include "AliEmcalJet.h"
42 #include "AliRhoParameter.h"
43 #include "AliLog.h"
44 #include "AliJetContainer.h"
45 #include "AliTrackContainer.h"
46 #include "AliAODTrack.h"
47 #include "AliPicoTrack.h"
48 #include "AliVParticle.h"
49 #include "TRandom3.h"
51 
53 
54 //________________________________________________________________________
56 {
57 // dummy destructor
58 }
59 
60 //________________________________________________________________________
62 {
63 // dummy destructor
64 }
65 
66 //________________________________________________________________________
68 {
69 // dummy destructor
70 }
71 
72 //________________________________________________________________________
74 {
75 // dummy destructor
76 }
77 
81 //________________________________________________________________________
84  fJetsCont(0),
85  fTracksCont(0),
86  fJetsTree(0),
87  fJetsTreeBuffer(0),
88  fRandom(0),
89  fCurrentNJetsInEvents(0),
90  fCurrentJetTypeHM(0),
91  fCurrentJetTypeIC(0),
92  fCurrentLeadingJet(0),
93  fCurrentSubleadingJet(0),
94  fCurrentInitialParton1(0),
95  fCurrentInitialParton2(0),
96  fCurrentInitialParton1Type(0),
97  fCurrentInitialParton2Type(0),
98  fFoundIC(kFALSE),
99  fExtractionCutMinCent(-1),
100  fExtractionCutMaxCent(-1),
101  fExtractionCutUseIC(kFALSE),
102  fExtractionCutUseHM(kFALSE),
103  fExtractionCutMinPt(0.),
104  fExtractionCutMaxPt(200.),
105  fExtractionPercentage(1.0),
106  fExtractionListPIDsHM(),
107  fHadronMatchingRadius(0.5),
108  fInitialCollisionMatchingRadius(0.3),
109  fTruthParticleArrayName("mcparticles"),
110  fSecondaryVertexMaxChi2(1e10),
111  fSecondaryVertexMaxDispersion(0.02),
112  fAddPIDSignal(kFALSE)
113 {
114  // Default constructor.
115  SetMakeGeneralHistograms(kTRUE);
116  fRandom = new TRandom3(0);
117 }
118 
119 //________________________________________________________________________
121  AliAnalysisTaskEmcalJet(name, kTRUE),
122  fJetsCont(0),
123  fTracksCont(0),
124  fJetsTree(0),
125  fJetsTreeBuffer(0),
126  fRandom(0),
127  fCurrentNJetsInEvents(0),
128  fCurrentJetTypeHM(0),
129  fCurrentJetTypeIC(0),
130  fCurrentLeadingJet(0),
131  fCurrentSubleadingJet(0),
132  fCurrentInitialParton1(0),
133  fCurrentInitialParton2(0),
134  fCurrentInitialParton1Type(0),
135  fCurrentInitialParton2Type(0),
136  fFoundIC(kFALSE),
137  fExtractionCutMinCent(-1),
138  fExtractionCutMaxCent(-1),
139  fExtractionCutUseIC(kFALSE),
140  fExtractionCutUseHM(kFALSE),
141  fExtractionCutMinPt(0.),
142  fExtractionCutMaxPt(200.),
143  fExtractionPercentage(1.0),
144  fExtractionListPIDsHM(),
145  fHadronMatchingRadius(0.5),
146  fInitialCollisionMatchingRadius(0.3),
147  fTruthParticleArrayName("mcparticles"),
148  fSecondaryVertexMaxChi2(1e10),
149  fSecondaryVertexMaxDispersion(0.02),
150  fAddPIDSignal(kFALSE)
151 {
152  // Default constructor.
154  fRandom = new TRandom3(0);
155 }
156 
157 //________________________________________________________________________
159 {
160  // Destructor.
161 }
162 
163 //________________________________________________________________________
165 {
167 
168  // ### Basic container settings
170  if(!fJetsCont)
171  AliFatal("Jet input container not found!");
172  fJetsCont->PrintCuts();
174  if(!fTracksCont)
175  AliFatal("Particle input container not found attached to jets!");
176 
177  //if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
178 
179  // ### Add control histograms (already some created in base task)
180  AddHistogram2D<TH2D>("hJetCount", "Number of jets in acceptance vs. centrality", "COLZ", 100, 0., 100., 100, 0, 100, "N Jets","Centrality", "dN^{Events}/dN^{Jets}");
181  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
182  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
183 
184  AddHistogram2D<TH2D>("hJetPtRaw", "Jets p_{T} distribution (raw)", "COLZ", 300, 0., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
185  AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution (background subtracted)", "COLZ", 400, -100., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
186  AddHistogram2D<TH2D>("hJetPhiEta", "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
187  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
188 
189  AddHistogram2D<TH2D>("hConstituentPt", "Jet constituent p_{T} distribution", "COLZ", 400, 0., 300., 100, 0, 100, "p_{T, const} (GeV/c)", "Centrality", "dN^{Const}/dp_{T}");
190  AddHistogram2D<TH2D>("hConstituentPhiEta", "Jet constituent relative #phi/#eta distribution", "COLZ", 120, -0.6, 0.6, 120, -0.6, 0.6, "#Delta#phi", "#Delta#eta", "dN^{Const}/d#phi d#eta");
191 
192  TH1* tmpHisto = AddHistogram1D<TH1D>("hJetAcceptance", "Accepted jets", "", 5, 0, 5, "stage","N^{jets}/cut");
193  tmpHisto->GetXaxis()->SetBinLabel(1, "Before cuts");
194  tmpHisto->GetXaxis()->SetBinLabel(2, "Centrality");
195  tmpHisto->GetXaxis()->SetBinLabel(3, "p_{T}");
196  tmpHisto->GetXaxis()->SetBinLabel(4, "PID");
197  tmpHisto->GetXaxis()->SetBinLabel(5, "Extraction percentage");
198 
199  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
200 }
201 
202 
203 //________________________________________________________________________
205 
207 
208  // ### Prepare the jet tree
209  fJetsTree = new TTree("ExtractedJets", "ExtractedJets");
210  fJetsTree->Branch("Jets", "AliBasicJet", &fJetsTreeBuffer, 1000);
211  fOutput->Add(fJetsTree);
212 }
213 
214 //________________________________________________________________________
216 {
217  // ######## Check if cuts for extraction have been passed
218 
219  // Extraction pt, jet PID (according to initial parton, according to HM), centrality
220  FillHistogram("hJetAcceptance", 0.5);
221 
222  // ### CENTRALITY
223  if(fExtractionCutMinCent != -1 && (fCent < fExtractionCutMinCent || fCent >= fExtractionCutMaxCent))
224  return kFALSE;
225  FillHistogram("hJetAcceptance", 1.5);
226 
227  // ### PT
228  if( ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) < fExtractionCutMinPt) || ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) >= fExtractionCutMaxPt) )
229  return kFALSE;
230  FillHistogram("hJetAcceptance", 2.5);
231 
232  Bool_t passedCutPID = kTRUE;
233  // ### PID (from initial collision)
235  {
236  passedCutPID = kFALSE;
237  for(Int_t i=0; i<fExtractionListPIDsIC.size(); i++)
238  {
240  passedCutPID = kTRUE;
241  }
242  }
243  // ### PID (from hadron matching)
244  else if(fExtractionCutUseHM)
245  {
246  passedCutPID = kFALSE;
247  for(Int_t i=0; i<fExtractionListPIDsHM.size(); i++)
248  {
250  passedCutPID = kTRUE;
251  }
252  }
253 
254  // Note: When either of the cuts is passed, accepted it
255  if(!passedCutPID)
256  return kFALSE;
257  FillHistogram("hJetAcceptance", 3.5);
258 
259  // Discard jets statistically
260  if(fRandom->Rndm() >= fExtractionPercentage)
261  return kFALSE;
262  FillHistogram("hJetAcceptance", 4.5);
263 
265  return kTRUE;
266 }
267 
268 //________________________________________________________________________
270 {
271  // Load vertex if possible
272  Double_t vtxX = 0;
273  Double_t vtxY = 0;
274  Double_t vtxZ = 0;
275  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
276  if(!myVertex && MCEvent())
277  myVertex = MCEvent()->GetPrimaryVertex();
278  if(myVertex)
279  {
280  vtxX = myVertex->GetX();
281  vtxY = myVertex->GetY();
282  vtxZ = myVertex->GetZ();
283  }
284 
285  // Get event ID from header
286  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
287  Long64_t eventID = 0;
288  if(eventIDHeader)
289  eventID = eventIDHeader->GetEventIdAsLong();
290 
291  // ### Actually add the basic jet
292  AliBasicJet basicJet(jet->Eta(), jet->Phi(), jet->Pt(), jet->Charge(), fJetsCont->GetJetRadius(), jet->Area(), fCurrentJetTypeIC, fCurrentJetTypeHM, fJetsCont->GetRhoVal(), InputEvent()->GetMagneticField(), vtxX, vtxY, vtxZ, eventID, fCent);
293 
294  // Add constituents
295  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
296  {
297  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
298  if(!particle) continue;
299 
300  // Secondary vertex analysis
301  Double_t z = GetTrackImpactParameter(myVertex, dynamic_cast<AliAODTrack*>(particle));
302 
303  basicJet.AddJetConstituent(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge(), particle->Xv(), particle->Yv(), particle->Zv(), z);
304  AddPIDInformation(particle, *basicJet.GetJetConstituent(basicJet.GetNumbersOfConstituents()-1));
305  }
306 
307  // Add secondary vertices to jet (from AOD.VertexingHF)
308  AddSecondaryVertices(myVertex, jet, basicJet);
309 
310  // Field currently not used
311  basicJet.SetTruePt(0);
312 
313  fJetsTreeBuffer = &basicJet;
314  fJetsTree->Fill();
315 }
316 
317 //________________________________________________________________________
319 {
320  // ### Event control plots
323  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
324 }
325 
326 //________________________________________________________________________
328 {
329  // ### Jet control plots
330  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
331  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
332  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
333  FillHistogram("hJetArea", jet->Area(), fCent);
334 
335  // ### Jet constituent plots
336  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
337  {
338  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
339  if(!jetConst) continue;
340 
341  // Constituent eta/phi (relative to jet)
342  Double_t deltaEta = jet->Eta() - jetConst->Eta();
343  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
344  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
345  deltaPhi = -deltaPhi;
346 
347  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
348  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
349  }
350 }
351 
352 //________________________________________________________________________
354 {
356 
357  // ### Jet loop
358  fJetsCont->ResetCurrentID();
359  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
360  {
362 
363  if(!IsJetSelected(jet))
364  continue;
365 
367  AddJetToTree(jet);
368  }
369 
371 
372  return kTRUE;
373 }
374 
375 //________________________________________________________________________
377 {
378  if(!vtx)
379  return;
380 
381  // Create ESD vertex from the existing AliVVertex
382  Double_t vtxPos[3] = {vtx->GetX(), vtx->GetY(), vtx->GetZ()};
383  Double_t covMatrix[6] = {0};
384  vtx->GetCovarianceMatrix(&covMatrix[0]);
385  AliESDVertex* esdVtx = new AliESDVertex(&vtxPos[0], &covMatrix[0], vtx->GetChi2(), vtx->GetNContributors());
386 
387  // Load HF vertex branch from AOD event, if possible
388  TClonesArray* secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
389  if(!secVertexArr)
390  return;
391 
392  // Loop over all potential secondary vertices
393  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
394  {
395  AliAODVertex* vtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
396  // Only select three-prong vertices
397  if((strcmp(vtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
398  continue;
399 
400  // Calculate vtx distance
401  Double_t effX = vtx->GetX() - esdVtx->GetX();
402  Double_t effY = vtx->GetY() - esdVtx->GetY();
403  Double_t effZ = vtx->GetZ() - esdVtx->GetZ();
404  Double_t vtxDistance = TMath::Sqrt(effX*effX + effY*effY + effZ*effZ);
405 
406  // Vertex dispersion
407  Double_t dispersion = 0.;
408  for(Int_t j=0; j<vtx->GetNDaughters(); j++)
409  {
410  if(!vtx->GetDaughter(j))
411  continue;
412  Double_t impact = GetTrackImpactParameter(vtx, dynamic_cast<AliAODTrack*>(vtx->GetDaughter(j)));
413  dispersion += impact*impact;
414  }
415  dispersion = TMath::Sqrt(dispersion);
416 
417  // Add secondary vertices if they fulfill the conditions
418  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(vtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
419  continue;
420 
421  basicJet.AddSecondaryVertex(vtx->GetX(), vtx->GetY(), vtx->GetZ(), TMath::Abs(vtx->GetChi2perNDF()), dispersion);
422  }
423 
424  delete esdVtx;
425 }
426 
427 //________________________________________________________________________
428 Double_t AliAnalysisTaskJetExtractor::GetTrackImpactParameter(const AliVVertex* vtx, AliAODTrack* track)
429 {
430  Double_t z = 0;
431  // Impact parameter z
432  if (track)
433  {
434  Double_t d0rphiz[2],covd0[3];
435  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),9999.,d0rphiz,covd0);
436  if(isDCA)
437  z = d0rphiz[1]; // impact parameter z
438  }
439 
440  return z;
441 }
442 
443 //________________________________________________________________________
445 {
447 }
448 
449 //________________________________________________________________________
451 {
455 }
456 
457 //________________________________________________________________________
459 {
460  typeIC = -1;
461  typeHM = -1;
462 
463  // Set type if initial parton information is available
464  if(fFoundIC)
465  {
466  typeIC = 0;
467  if(jet==fCurrentInitialParton1)
469  else if (jet==fCurrentInitialParton2)
471  }
472 
473  // Do hadron matching jet type tagging using mcparticles
474  // ... if not explicitly deactivated
475  if (fTruthParticleArrayName != "")
476  {
477  TClonesArray* fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
478  if(!fTruthParticleArray)
479  return;
480 
481  typeHM = 0;
482  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
483  {
484  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
485  if(!part) continue;
486 
487  // Check if particle is in a radius around the jet
488  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
490  continue;
491 
492  // Check if the particle has beauty, charm or strangeness
493  // If it has beauty, we are done (exclusive definition)
494  Int_t absPDG = TMath::Abs(part->PdgCode());
495  // Particle has beauty
496  if ((absPDG > 500 && absPDG < 600) || (absPDG > 5000 && absPDG < 6000))
497  {
498  typeHM = 5; // beauty
499  break;
500  }
501  // Particle has charm
502  else if ((absPDG > 400 && absPDG < 500) || (absPDG > 4000 && absPDG < 5000))
503  typeHM = 4; // charm
504  // Particle has strangeness: Only search for strangeness, if charm was not already found
505  else if (typeHM != 4 && (absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
506  typeHM = 3; // strange
507  }
508  }
509  // As fallback, the MC stack will be tried
510  else if(MCEvent() && (MCEvent()->Stack()))
511  {
512  typeHM = 0;
513  AliStack* stack = MCEvent()->Stack();
514  // Go through the whole particle stack
515  for(Int_t i=0; i<stack->GetNtrack(); i++)
516  {
517  TParticle *part = stack->Particle(i);
518  if(!part) continue;
519 
520  // Check if particle is in a radius around the jet
521  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
523  continue;
524 
525  // Check if the particle has beauty, charm or strangeness
526  // If it has beauty, we are done (exclusive definition)
527  Int_t absPDG = TMath::Abs(part->GetPdgCode());
528  // Particle has beauty
529  if ((absPDG > 500 && absPDG < 600) || (absPDG > 5000 && absPDG < 6000))
530  {
531  typeHM = 5; // beauty
532  break;
533  }
534  // Particle has charm
535  else if ((absPDG > 400 && absPDG < 500) || (absPDG > 4000 && absPDG < 5000))
536  typeHM = 4; // charm
537  // Particle has strangeness: Only search for strangeness, if charm was not already found
538  else if (typeHM != 4 && (absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
539  typeHM = 3; // strange
540  }
541  }
542 }
543 
544 //________________________________________________________________________
546 {
547  // On demand and if we have AODs, add the particle PID signal object
548  // Add the truth values if available
549  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
550 
551  if(!fAddPIDSignal)
552  return;
553  if(!aodtrack)
554  return;
555 
556  // Get AOD value from reco
557  Int_t recoPID = aodtrack->GetMostProbablePID();
558  Int_t truthPID = 9;
559 
560  // Get truth values if we are on MC
561  TClonesArray* fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
562  if(fTruthParticleArray)
563  {
564  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
565  {
566  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
567  if(!mcParticle) continue;
568 
569  if (mcParticle->GetLabel() == particle->GetLabel())
570  {
571  // Use same convention as PID in AODs
572  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
573  truthPID = 4;
574  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
575  truthPID = 2;
576  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
577  truthPID = 3;
578  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
579  truthPID = 0;
580  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
581  truthPID = 1;
582  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
583  truthPID = 5;
584  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
585  truthPID = 6;
586  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
587  truthPID = 7;
588  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
589  truthPID = 8;
590  else
591  truthPID = 9;
592 
593  break;
594  }
595  }
596  }
597 
598  AliAODPid* pidObj = aodtrack->GetDetPid();
599  constituent.SetPIDSignal(pidObj->GetITSsignal(), pidObj->GetTPCsignal(), pidObj->GetTOFsignal(), pidObj->GetTRDsignal(), truthPID, recoPID);
600 }
601 
602 //________________________________________________________________________
603 void AliAnalysisTaskJetExtractor::GetLeadingJets(const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
604 {
605  // Customized from AliJetContainer::GetLeadingJet()
606  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
607 
608  TString option(opt);
609  option.ToLower();
610 
611  jetLeading = 0;
612  jetSubLeading = 0;
613 
614  fJetsCont->ResetCurrentID();
615  Double_t tmpLeadingPt = 0;
616  Double_t tmpSubleadingPt = 0;
617 
618  if (option.Contains("rho")) {
619  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
620  if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpLeadingPt )
621  {
622  jetSubLeading = jetLeading;
623  jetLeading = jet;
624  tmpSubleadingPt = tmpLeadingPt;
625  tmpLeadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
626  }
627  else if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpSubleadingPt )
628  {
629  jetSubLeading = jet;
630  tmpSubleadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
631  }
632  }
633  }
634  else {
635  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
636  if ( (jet->Pt()) > tmpLeadingPt )
637  {
638  jetSubLeading = jetLeading;
639  jetLeading = jet;
640  tmpSubleadingPt = tmpLeadingPt;
641  tmpLeadingPt = jet->Pt();
642  }
643  else if ( (jet->Pt()) > tmpSubleadingPt )
644  {
645  jetSubLeading = jet;
646  tmpSubleadingPt = jet->Pt();
647  }
648  }
649  }
650 }
651 
652 //________________________________________________________________________
654 {
655  // Get the initial parton infromation
656  Double_t initialParton1_eta = -999.;
657  Double_t initialParton1_phi = -999.;
658  Int_t initialParton1_pdg = 0;
659 
660  Double_t initialParton2_eta = -999.;
661  Double_t initialParton2_phi = -999.;
662  Int_t initialParton2_pdg = 0;
663 
664  // If available, use stack as input
665  if(MCEvent() && (MCEvent()->Stack()))
666  {
667  AliStack* stack = MCEvent()->Stack();
668  TParticle* parton1 = stack->Particle(6);
669  TParticle* parton2 = stack->Particle(7);
670  if(parton1)
671  {
672  initialParton1_eta = parton1->Eta();
673  initialParton1_phi = parton1->Phi();
674  initialParton1_pdg = TMath::Abs(parton1->GetPdgCode());
675  }
676  if(parton2)
677  {
678  initialParton2_eta = parton2->Eta();
679  initialParton2_phi = parton2->Phi();
680  initialParton2_pdg = TMath::Abs(parton2->GetPdgCode());
681  }
682  }
683  // Otherwise, PYTHIA object information
684  else if (fPythiaInfo)
685  {
686  initialParton1_eta = fPythiaInfo->GetPartonEta6();
687  initialParton1_phi = fPythiaInfo->GetPartonPhi6();
688  initialParton1_pdg = fPythiaInfo->GetPartonFlag6();
689  initialParton2_eta = fPythiaInfo->GetPartonEta7();
690  initialParton2_phi = fPythiaInfo->GetPartonPhi7();
691  initialParton2_pdg = fPythiaInfo->GetPartonFlag7();
692  }
693  // or direct particle level information from mcparticles branch
694  else if (TClonesArray* fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data())))
695  {
696  AliAODMCParticle* parton1 = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(6));
697  AliAODMCParticle* parton2 = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(7));
698  if(parton1)
699  {
700  initialParton1_eta = parton1->Eta();
701  initialParton1_phi = parton1->Phi();
702  initialParton1_pdg = TMath::Abs(parton1->GetPdgCode());
703  }
704  if(parton2)
705  {
706  initialParton2_eta = parton2->Eta();
707  initialParton2_phi = parton2->Phi();
708  initialParton2_pdg = TMath::Abs(parton2->GetPdgCode());
709  }
710  }
711  // No initial collision partons found, return
712  else
713  {
714  fFoundIC = kFALSE;
719  return;
720  }
721 
722  fFoundIC = kTRUE;
723 
724  Double_t bestMatchDeltaR1 = 999.;
725  Double_t bestMatchDeltaR2 = 999.;
726 
727  // #### Find via geometrical matching the jet connected to the initial collision
728  fJetsCont->ResetCurrentID();
729  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
730  {
731  Double_t deltaEta1 = TMath::Abs(jet->Eta()-initialParton1_eta);
732  Double_t deltaEta2 = TMath::Abs(jet->Eta()-initialParton2_eta);
733  Double_t deltaPhi1 = TMath::Min(TMath::Abs(jet->Phi()-initialParton1_phi),TMath::TwoPi() - TMath::Abs(jet->Phi()-initialParton1_phi));
734  Double_t deltaPhi2 = TMath::Min(TMath::Abs(jet->Phi()-initialParton2_phi),TMath::TwoPi() - TMath::Abs(jet->Phi()-initialParton2_phi));
735 
736  Double_t deltaR1 = TMath::Sqrt(deltaEta1*deltaEta1 + deltaPhi1*deltaPhi1);
737  Double_t deltaR2 = TMath::Sqrt(deltaEta2*deltaEta2 + deltaPhi2*deltaPhi2);
738 
739  if(deltaR1 < bestMatchDeltaR1)
740  {
741  bestMatchDeltaR1 = deltaR1;
743  fCurrentInitialParton1Type = initialParton1_pdg;
744  }
745  if(deltaR2 < bestMatchDeltaR2)
746  {
747  bestMatchDeltaR2 = deltaR2;
749  fCurrentInitialParton2Type = initialParton2_pdg;
750  }
751  }
752 
753  // Discard matched jet that are to far away
754  if(bestMatchDeltaR1 > fInitialCollisionMatchingRadius)
756  if(bestMatchDeltaR2 > fInitialCollisionMatchingRadius)
758 }
759 
760 
761 
762 //########################################################################
763 // HELPERS
764 //########################################################################
765 
766 
767 //________________________________________________________________________
768 inline void AliAnalysisTaskJetExtractor::FillHistogram(const char * key, Double_t x)
769 {
770  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
771  if(!tmpHist)
772  {
773  AliError(Form("Cannot find histogram <%s> ",key)) ;
774  return;
775  }
776 
777  tmpHist->Fill(x);
778 }
779 
780 //________________________________________________________________________
782 {
783  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
784  if(!tmpHist)
785  {
786  AliError(Form("Cannot find histogram <%s> ",key));
787  return;
788  }
789 
790  if (tmpHist->IsA()->GetBaseClass("TH1"))
791  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
792  else if (tmpHist->IsA()->GetBaseClass("TH2"))
793  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
794 }
795 
796 //________________________________________________________________________
798 {
799  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
800  if(!tmpHist)
801  {
802  AliError(Form("Cannot find histogram <%s> ",key));
803  return;
804  }
805 
806  tmpHist->Fill(x,y,add);
807 }
808 
809 //________________________________________________________________________
811 {
812  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
813  if(!tmpHist)
814  {
815  AliError(Form("Cannot find histogram <%s> ",key));
816  return;
817  }
818 
819  if(add)
820  tmpHist->Fill(x,y,z,add);
821  else
822  tmpHist->Fill(x,y,z);
823 }
824 
825 
826 //________________________________________________________________________
827 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
828 {
829  T* tmpHist = new T(name, title, xBins, xMin, xMax);
830 
831  tmpHist->GetXaxis()->SetTitle(xTitle);
832  tmpHist->GetYaxis()->SetTitle(yTitle);
833  tmpHist->SetOption(options);
834  tmpHist->SetMarkerStyle(kFullCircle);
835  tmpHist->Sumw2();
836 
837  fOutput->Add(tmpHist);
838 
839  return tmpHist;
840 }
841 
842 //________________________________________________________________________
843 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
844 {
845  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
846  tmpHist->GetXaxis()->SetTitle(xTitle);
847  tmpHist->GetYaxis()->SetTitle(yTitle);
848  tmpHist->GetZaxis()->SetTitle(zTitle);
849  tmpHist->SetOption(options);
850  tmpHist->SetMarkerStyle(kFullCircle);
851  tmpHist->Sumw2();
852 
853  fOutput->Add(tmpHist);
854 
855  return tmpHist;
856 }
857 
858 //________________________________________________________________________
859 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram3D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, Int_t zBins, Double_t zMin, Double_t zMax, const char* xTitle, const char* yTitle, const char* zTitle)
860 {
861  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
862  tmpHist->GetXaxis()->SetTitle(xTitle);
863  tmpHist->GetYaxis()->SetTitle(yTitle);
864  tmpHist->GetZaxis()->SetTitle(zTitle);
865  tmpHist->SetOption(options);
866  tmpHist->SetMarkerStyle(kFullCircle);
867  tmpHist->Sumw2();
868 
869  fOutput->Add(tmpHist);
870 
871  return tmpHist;
872 }
873 
874 //________________________________________________________________________
876 {
877  // Called once at the end of the analysis.
878 }
879 
Int_t fExtractionCutMinCent
Extraction cut: minimum centrality.
void CalculateJetType(AliEmcalJet *jet, Int_t &typeIC, Int_t &typeHM)
void AddPIDInformation(AliVParticle *particle, AliBasicJetConstituent &constituent)
AliEmcalJet * fCurrentInitialParton1
! jet that matches the initial parton 1 (PYTHIA)
Short_t Charge() const
Definition: AliEmcalJet.h:110
Double_t Area() const
Definition: AliEmcalJet.h:117
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
TString fTruthParticleArrayName
Array name of MC particles in event (mcparticles)
const char * title
Definition: MakeQAPdf.C:26
Float_t GetPartonEta6() const
AliEmcalPythiaInfo * fPythiaInfo
!event parton info
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
Float_t GetPartonEta7() const
Double_t fExtractionCutMinPt
Extraction cut: minimum pT.
Double_t Eta() const
Definition: AliEmcalJet.h:108
long long Long64_t
Definition: External.C:43
T * AddHistogram3D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, Int_t zBins=100, Double_t zMin=0.0, Double_t zMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
AliEmcalJet * fCurrentLeadingJet
! leading jet (calculated event-by-event)
Double_t Phi() const
Definition: AliEmcalJet.h:104
Int_t fExtractionCutMaxCent
Extraction cut: maximum centrality.
Container with name, TClonesArray and cuts for particles.
Double_t fSecondaryVertexMaxChi2
Max chi2 of secondary vertex (others will be discarded)
Float_t GetPartonPhi7() const
Double_t fSecondaryVertexMaxDispersion
Max dispersion of secondary vertex (others will be discarded)
AliTrackContainer * fTracksCont
! Tracks
void AddSecondaryVertices(const AliVVertex *vtx, AliEmcalJet *jet, AliBasicJet &basicJet)
Simple class containing basic information for a constituent.
T * AddHistogram2D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
void AddSecondaryVertex(Float_t vx, Float_t vy, Float_t vz, Float_t chi2=0, Float_t dispersion=0)
Int_t fCurrentInitialParton2Type
type of initial parton 2
Int_t GetPartonFlag7() const
Simple class containing basic information for a jet.
Bool_t fExtractionCutUseIC
Extraction cut: Use IC PID.
Float_t GetPartonPhi6() const
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
AliEmcalJet * fCurrentInitialParton2
! jet that matches the initial parton 2 (PYTHIA)
void FillHistogram3D(const char *key, Double_t x, Double_t y, Double_t z, Double_t add=0)
int Int_t
Definition: External.C:63
Int_t fCurrentNJetsInEvents
number of jets in event
void * fJetsTreeBuffer
! buffer for one jet (that will be saved to the tree)
TRandom3 * fRandom
random number generator
T * AddHistogram1D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis")
Bool_t fFoundIC
status var showing that IC has been found
Double_t fCent
!event centrality
Double_t fExtractionPercentage
Percentage of extracted jets.
std::vector< Int_t > fExtractionListPIDsIC
list of PIDs (initial collision) that will be accepted
AliEmcalJet * fCurrentSubleadingJet
! subleading jet (calculated event-by-event)
Double_t GetTrackImpactParameter(const AliVVertex *vtx, AliAODTrack *track)
AliEmcalJet * GetNextAcceptJet()
Double_t Pt() const
Definition: AliEmcalJet.h:96
void FillHistogram(const char *key, Double_t x)
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
Definition: External.C:220
Bool_t fExtractionCutUseHM
Extraction cut: Use hadron matching PID.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
std::vector< Int_t > fExtractionListPIDsHM
list of PIDs (hadron matching) that will be accepted
Bool_t fAddPIDSignal
Add pid signal to each jet constituent.
Int_t fCurrentInitialParton1Type
type of initial parton 1
Double_t fExtractionCutMaxPt
Extraction cut: minimum pT.
void SetMakeGeneralHistograms(Bool_t g)
Double_t fHadronMatchingRadius
Matching radius to search for beauty/charm hadrons around jet.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Declaration of class AliEmcalPythiaInfo.
Int_t fCurrentJetTypeHM
current jet type from hadron matching
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void GetLeadingJets(const char *opt, AliEmcalJet *&jetLeading, AliEmcalJet *&jetSubLeading)
Int_t GetNAcceptedParticles() const
void SetPIDSignal(Float_t its, Float_t tpc, Float_t tof, Float_t trd, Short_t truthPID, Short_t recoPID)
TTree * fJetsTree
! Jets that will be saved to a tree (optionally)
Int_t fCurrentJetTypeIC
current jet type from initial collision (PYTHIA)
void FillJetControlHistograms(AliEmcalJet *jet)
Double_t fInitialCollisionMatchingRadius
Matching radius to find a jet of the IC.
Definition: External.C:196