AliPhysics  84d6c4a (84d6c4a)
AliAnalysisTaskSoftDropResponse.cxx
Go to the documentation of this file.
1 /*************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
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 // Soft Drop observables response taks based on the AliJetResponseMaker task
17 //
18 // Author: Kirill Lapidus, Yale University, kirill.lapdidus@cern.ch
19 
20 #include <THnSparse.h>
21 
23 #include "AliAnalysisManager.h"
26 
28 
29 //________________________________________________________________________
32  fZ2gAxis(0)
33 {
34  // Default constructor.
35 
36  SetMakeGeneralHistograms(kTRUE);
37 }
38 
39 //________________________________________________________________________
41  AliJetResponseMaker(name),
42  fZ2gAxis(0)
43 {
44  // Standard constructor.
45 
47 }
48 
49 
50 //________________________________________________________________________
52 {
53  // Allocate THnSparse histograms.
54 
55  TString title[25]= {""};
56  Int_t nbins[25] = {0};
57  Double_t min[25] = {0.};
58  Double_t max[25] = {0.};
59  Int_t dim = 0;
60 
61  title[dim] = "#phi";
62  nbins[dim] = fNbins/4;
63  min[dim] = 0;
64  max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
65  dim++;
66 
67  title[dim] = "#eta";
68  nbins[dim] = fNbins/4;
69  min[dim] = -1;
70  max[dim] = 1;
71  dim++;
72 
73  title[dim] = "p_{T}";
74  nbins[dim] = fNbins;
75  min[dim] = 0;
76  max[dim] = 250;
77  dim++;
78 
79  title[dim] = "A_{jet}";
80  nbins[dim] = fNbins/4;
81  min[dim] = 0;
82  max[dim] = 1.5;
83  dim++;
84 
85  Int_t dim1 = dim, dim2 = dim;
86 
87  if (fIsJet1Rho) {
88  title[dim1] = "p_{T}^{corr}";
89  nbins[dim1] = fNbins*2;
90  min[dim1] = -250;
91  max[dim1] = 250;
92  dim1++;
93  }
94 
95  if (fIsEmbedded) {
96  title[dim1] = "p_{T}^{MC}";
97  nbins[dim1] = fNbins;
98  min[dim1] = 0;
99  max[dim1] = 250;
100  dim1++;
101  }
102 
103  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
104  for (Int_t i = 0; i < dim1; i++) fHistJets1->GetAxis(i)->SetTitle(title[i]);
105  fOutput->Add(fHistJets1);
106 
107  if (fIsJet2Rho) {
108  title[dim2] = "p_{T}^{corr}";
109  nbins[dim2] = fNbins*2;
110  min[dim2] = -250;
111  max[dim2] = 250;
112  dim2++;
113  }
114 
115  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
116  for (Int_t i = 0; i < dim2; i++) fHistJets2->GetAxis(i)->SetTitle(title[i]);
117  fOutput->Add(fHistJets2);
118 
119  // Matching
120 
121  dim = 0;
122 
123  title[dim] = "p_{T,1}";
124  nbins[dim] = fNbins;
125  min[dim] = 0;
126  max[dim] = 250;
127  dim++;
128 
129  title[dim] = "p_{T,2}";
130  nbins[dim] = fNbins;
131  min[dim] = 0;
132  max[dim] = 250;
133  dim++;
134 
135  title[dim] = "A_{jet,1}";
136  nbins[dim] = fNbins/4;
137  min[dim] = 0;
138  max[dim] = 1.5;
139  dim++;
140 
141  title[dim] = "A_{jet,2}";
142  nbins[dim] = fNbins/4;
143  min[dim] = 0;
144  max[dim] = 1.5;
145  dim++;
146 
147  title[dim] = "distance";
148  nbins[dim] = fNbins/2;
149  min[dim] = 0;
150  max[dim] = 1.2;
151  dim++;
152 
153  if (fDeltaPtAxis) {
154  title[dim] = "#deltaA_{jet}";
155  nbins[dim] = fNbins/2;
156  min[dim] = -1;
157  max[dim] = 1;
158  dim++;
159 
160  title[dim] = "#deltap_{T}";
161  nbins[dim] = fNbins*2;
162  min[dim] = -250;
163  max[dim] = 250;
164  dim++;
165  }
166  if (fIsJet1Rho) {
167  title[dim] = "p_{T,1}^{corr}";
168  nbins[dim] = fNbins*2;
169  min[dim] = -250;
170  max[dim] = 250;
171  dim++;
172  }
173  if (fIsJet2Rho) {
174  title[dim] = "p_{T,2}^{corr}";
175  nbins[dim] = fNbins*2;
176  min[dim] = -250;
177  max[dim] = 250;
178  dim++;
179  }
180  if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
181  title[dim] = "#deltap_{T}^{corr}";
182  nbins[dim] = fNbins*2;
183  min[dim] = -250;
184  max[dim] = 250;
185  dim++;
186  }
187  if (fDeltaEtaDeltaPhiAxis) {
188  title[dim] = "#delta#eta";
189  nbins[dim] = fNbins/2;
190  min[dim] = -1;
191  max[dim] = 1;
192  dim++;
193 
194  title[dim] = "#delta#phi";
195  nbins[dim] = fNbins/2;
196  min[dim] = -TMath::Pi()/2;
197  max[dim] = TMath::Pi()*3/2;
198  dim++;
199  }
200  if (fIsEmbedded) {
201  title[dim] = "p_{T,1}^{MC}";
202  nbins[dim] = fNbins;
203  min[dim] = 0;
204  max[dim] = 250;
205  dim++;
206 
207  if (fDeltaPtAxis) {
208  title[dim] = "#deltap_{T}^{MC}";
209  nbins[dim] = fNbins*2;
210  min[dim] = -250;
211  max[dim] = 250;
212  dim++;
213  }
214  }
215 
216  if (fZgAxis) {
217  title[dim] = "Z_{g,1}";
218  nbins[dim] = 20;
219  min[dim] = 0.0;
220  max[dim] = 1.0;
221  dim++;
222  title[dim] = "Z_{g,2}";
223  nbins[dim] = 20;
224  min[dim] = 0.0;
225  max[dim] = 1.0;
226  dim++;
227  }
228 
229  if (fZ2gAxis) {
230  title[dim] = "Z2_{g,1}";
231  nbins[dim] = 20;
232  min[dim] = 0.0;
233  max[dim] = 1.0;
234  dim++;
235  title[dim] = "Z2_{g,2}";
236  nbins[dim] = 20;
237  min[dim] = 0.0;
238  max[dim] = 1.0;
239  dim++;
240  }
241 
242  if (fdRAxis) {
243  title[dim] = "R_{g,1}";
244  nbins[dim] = 20;
245  min[dim] = 0.0;
246  max[dim] = 0.5;
247  dim++;
248  title[dim] = "R_{g,2}";
249  nbins[dim] = 20;
250  min[dim] = 0.0;
251  max[dim] = 0.5;
252  dim++;
253  }
254 
255  if (fPtgAxis) {
256  title[dim] = "p_{T,g,1}";
257  nbins[dim] = 16;
258  min[dim] = 0.0;
259  max[dim] = 160.0;
260  dim++;
261  title[dim] = "p_{T,g,2}";
262  nbins[dim] = 16;
263  min[dim] = 0.0;
264  max[dim] = 160.0;
265  dim++;
266  }
267 
268  if (fDBCAxis) {
269  title[dim] = "DBC_{1}";
270  nbins[dim] = 20;
271  min[dim] = 0.0;
272  max[dim] = 20.0;
273  dim++;
274  title[dim] = "DBC_{2}";
275  nbins[dim] = 20;
276  min[dim] = 0.0;
277  max[dim] = 20.0;
278  dim++;
279  }
280 
281  title[dim] = "nsdsteps_{jet,1}";
282  nbins[dim] = 10;
283  min[dim] = 1;
284  max[dim] = 10;
285  dim++;
286 
287  title[dim] = "nsdsteps_{jet,2}";
288  nbins[dim] = 10;
289  min[dim] = 1;
290  max[dim] = 10;
291  dim++;
292 
293  title[dim] = "nhsplits_{jet,1}";
294  nbins[dim] = 10;
295  min[dim] = 1;
296  max[dim] = 10;
297  dim++;
298 
299  title[dim] = "nhsplits_{jet,2}";
300  nbins[dim] = 10;
301  min[dim] = 1;
302  max[dim] = 10;
303  dim++;
304 
305  title[dim] = "Z^{2}_{g,1}";
306  nbins[dim] = 20;
307  min[dim] = 0.0;
308  max[dim] = 0.5;
309  dim++;
310 
311  title[dim] = "Z^{2}_{g,2}";
312  nbins[dim] = 20;
313  min[dim] = 0.0;
314  max[dim] = 0.5;
315  dim++;
316 
317  title[dim] = "R^{2}_{g,1}";
318  nbins[dim] = 20;
319  min[dim] = 0.0;
320  max[dim] = 0.5;
321  dim++;
322 
323  title[dim] = "R^{2}_{g,2}";
324  nbins[dim] = 20;
325  min[dim] = 0.0;
326  max[dim] = 0.5;
327  dim++;
328 
329  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
330 
331  for (Int_t i = 0; i < dim; i++) fHistMatching->GetAxis(i)->SetTitle(title[i]);
332 
333  fOutput->Add(fHistMatching);
334 }
335 
336 //________________________________________________________________________
338 {
339  // Fill histograms.
340 
341  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
342  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
343 
344  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
345 
346  AliEmcalJet* jet1 = 0;
347  AliEmcalJet* jet2 = 0;
348 
349  jets2->ResetCurrentID();
350  while ((jet2 = jets2->GetNextJet())) {
351 
352  AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
353 
354  if (jet2->Pt() < jets2->GetJetPtCut()) continue;
355 
356  UInt_t rejectionReason = 0;
357  if (jets2->AcceptJet(jet2, rejectionReason))
358  FillJetHisto(jet2, 2);
359  else
360  fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(rejectionReason), jet2->Pt());
361 
362  jet1 = jet2->MatchedJet();
363 
364  if (!jet1) continue;
365  rejectionReason = 0;
366  if (!jets1->AcceptJet(jet1, rejectionReason)) continue;
367  if (jet1->MCPt() < fMinJetMCPt) continue;
368 
369  Double_t d=-1, ce1=-1, ce2=-1;
370  if (jet2->GetMatchingType() == kGeometrical) {
371  if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
372  GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
373  else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
374  GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
375 
376  d = jet2->ClosestJetDistance();
377  }
378  else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
379  GetGeometricalMatchingLevel(jet1, jet2, d);
380 
381  ce1 = jet1->ClosestJetDistance();
382  ce2 = jet2->ClosestJetDistance();
383  }
384 
385  FillMatchingHistos(jet1, jet2, d, ce1, ce2);
386  }
387 
388  jets1->ResetCurrentID();
389  while ((jet1 = jets1->GetNextJet())) {
390  UInt_t rejectionReason = 0;
391  if (!jets1->AcceptJet(jet1, rejectionReason)) {
392  fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(rejectionReason), jet1->Pt());
393  continue;
394  }
395  if (jet1->MCPt() < fMinJetMCPt) continue;
396  AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
397 
398  FillJetHisto(jet1, 1);
399  }
400  return kTRUE;
401 }
402 
403 //________________________________________________________________________
405 {
406  AliJetContainer* jets1 = GetJetContainer(0);
407  AliJetContainer* jets2 = GetJetContainer(1);
408 
409  Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
410  Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
411 
412  //Float_t z2g1 = -1.0;
413  //Float_t z2g2 = -1.0;
414  //CalculateZg(jet1, 0.5, 1.5, z2g1);
415  //CalculateZg(jet2, 0.5, 1.5, z2g2);
416 
417  FillZgRgVectors(jet1);
418 
419  Int_t nsdsteps_jet1 = 0;
420  Int_t nhsplits_jet1 = fZg_values.size();
421  Float_t zg2_1 = 0.0;
422  Float_t rg2_1 = 0.0;
423  if (nhsplits_jet1 > 0) nsdsteps_jet1 = fSDsteps_values[0];
424  if (nhsplits_jet1 > 1) {
425  zg2_1 = fZg_values[1];
426  rg2_1 = fRg_values[1];
427  }
428 
429  FillZgRgVectors(jet2);
430 
431  Int_t nsdsteps_jet2 = 0;
432  Int_t nhsplits_jet2 = fZg_values.size();
433  Float_t zg2_2 = 0.0;
434  Float_t rg2_2 = 0.0;
435  if (nhsplits_jet2 > 0) nsdsteps_jet2 = fSDsteps_values[0];
436  if (nhsplits_jet2 > 1) {
437  zg2_2 = fZg_values[1];
438  rg2_2 = fRg_values[1];
439  }
440 
441  if (fHistoType==1) {
442  Double_t contents[25]={0};
443 
444  for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
445  TString title(fHistMatching->GetAxis(i)->GetTitle());
446  if (title=="p_{T,1}")
447  contents[i] = jet1->Pt();
448  else if (title=="p_{T,2}")
449  contents[i] = jet2->Pt();
450  else if (title=="A_{jet,1}")
451  contents[i] = jet1->Area();
452  else if (title=="A_{jet,2}")
453  contents[i] = jet2->Area();
454  else if (title=="distance")
455  contents[i] = d;
456  else if (title=="#deltaA_{jet}")
457  contents[i] = jet1->Area()-jet2->Area();
458  else if (title=="#deltap_{T}")
459  contents[i] = jet1->Pt()-jet2->Pt();
460  else if (title=="#delta#eta")
461  contents[i] = jet1->Eta()-jet2->Eta();
462  else if (title=="#delta#phi")
463  contents[i] = jet1->Phi()-jet2->Phi();
464  else if (title=="p_{T,1}^{corr}")
465  contents[i] = corrpt1;
466  else if (title=="p_{T,2}^{corr}")
467  contents[i] = corrpt2;
468  else if (title=="#deltap_{T}^{corr}")
469  contents[i] = corrpt1-corrpt2;
470  else if (title=="p_{T,1}^{MC}")
471  contents[i] = jet1->MCPt();
472  else if (title=="#deltap_{T}^{MC}")
473  contents[i] = jet1->MCPt()-jet2->Pt();
474  else if (title=="Z_{g,1}")
475  contents[i] = jet1->GetShapeProperties()->GetSoftDropZg();
476  else if (title=="Z_{g,2}")
477  contents[i] = jet2->GetShapeProperties()->GetSoftDropZg();
478  else if (title=="Z2_{g,1}")
479  contents[i] = 0.0;
480  else if (title=="Z2_{g,2}")
481  contents[i] = 0.0;
482  else if (title=="R_{g,1}")
483  contents[i] = jet1->GetShapeProperties()->GetSoftDropdR();
484  else if (title=="R_{g,2}")
485  contents[i] = jet2->GetShapeProperties()->GetSoftDropdR();
486  else if (title=="p_{T,g,1}")
487  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropPtfrac() )*( jet1->Pt() );
488  else if (title=="p_{T,g,2}")
489  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropPtfrac() )*( jet2->Pt() );
490  else if (title=="DBC_{1}")
491  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropDropCount() );
492  else if (title=="DBC_{2}")
493  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropDropCount() );
494  else if (title=="nsdsteps_{jet,1}")
495  contents[i] = nsdsteps_jet1;
496  else if (title=="nsdsteps_{jet,2}")
497  contents[i] = nsdsteps_jet2;
498  else if (title=="nhsplits_{jet,1}")
499  contents[i] = nhsplits_jet1;
500  else if (title=="nhsplits_{jet,2}")
501  contents[i] = nhsplits_jet2;
502  else if (title=="Z^{2}_{g,1}")
503  contents[i] = zg2_1;
504  else if (title=="Z^{2}_{g,2}")
505  contents[i] = zg2_2;
506  else if (title=="R^{2}_{g,1}")
507  contents[i] = rg2_1;
508  else if (title=="R^{2}_{g,2}")
509  contents[i] = rg2_2;
510  else
511  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
512  }
513 
514  fHistMatching->Fill(contents);
515  }
516 
517 }
518 
523  const char *ntracks1,
524  const char *nclusters1,
525  const char *njets1,
526  const char *nrho1,
527  const Double_t jetradius1,
528  const char *ntracks2,
529  const char *nclusters2,
530  const char *njets2,
531  const char *nrho2,
532  const Double_t jetradius2,
533  const Double_t jetptcut,
534  const Double_t jetareacut,
535  const Double_t jetBias,
536  const Int_t biasType,
538  const Double_t maxDistance1,
539  const Double_t maxDistance2,
540  const char *cutType,
541  const Int_t ptHardBin,
542  const Double_t minCent,
543  const Double_t maxCent,
544  const char *taskname,
545  const Bool_t biggerMatrix,
547  const Double_t maxTrackPt)
548 {
549  // Get the pointer to the existing analysis manager via the static access method.
550  //==============================================================================
551  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
552  if (!mgr)
553  {
554  ::Error("AddTaskJetResponseMaker", "No analysis manager to connect to.");
555  return NULL;
556  }
557 
558  // Check the analysis type using the event handlers connected to the analysis manager.
559  //==============================================================================
560  if (!mgr->GetInputEventHandler())
561  {
562  ::Error("AddTaskJetResponseMaker", "This task requires an input event handler");
563  return NULL;
564  }
565 
566  //-------------------------------------------------------
567  // Init the task and do settings
568  //-------------------------------------------------------
569 
570  TString name(Form("%s_%s_%s_Bias%d_BiasType%d_%s",taskname,njets1,njets2,(Int_t)floor(jetBias),biasType,cutType));
571 
572  if (minCent != -999 && maxCent != -999)
573  name += Form("_Cent%d_%d", (Int_t)floor(minCent), (Int_t)floor(maxCent));
574 
575  if (ptHardBin != -999)
576  name += Form("_PtHard%d", ptHardBin);
577 
578  AliAnalysisTaskSoftDropResponse* jetTask = address;
579  if (jetTask) {
580  new (jetTask) AliAnalysisTaskSoftDropResponse(name);
581  }
582  else {
583  jetTask = new AliAnalysisTaskSoftDropResponse(name);
584  }
585 
586  AliParticleContainer *trackCont1 = jetTask->AddParticleContainer(ntracks1);
587  AliClusterContainer *clusCont1 = jetTask->AddClusterContainer(nclusters1);
588  AliJetContainer *jetCont1 = jetTask->AddJetContainer(njets1, cutType, jetradius1);
589  jetCont1->SetRhoName(nrho1);
590  jetCont1->SetLeadingHadronType(biasType);
591  jetCont1->SetPtBiasJetTrack(jetBias);
592  jetCont1->SetPtBiasJetClus(jetBias);
593  jetCont1->SetJetPtCut(jetptcut);
594  jetCont1->SetPercAreaCut(jetareacut);
595  jetCont1->SetIsParticleLevel(kFALSE);
596  jetCont1->ConnectParticleContainer(trackCont1);
597  jetCont1->ConnectClusterContainer(clusCont1);
598  jetCont1->SetMaxTrackPt(maxTrackPt);
599 
600  AliParticleContainer *trackCont2 = jetTask->AddParticleContainer(ntracks2);
601  trackCont2->SetParticlePtCut(0);
602  AliClusterContainer *clusCont2 = jetTask->AddClusterContainer(nclusters2);
603  AliJetContainer *jetCont2 = jetTask->AddJetContainer(njets2, cutType, jetradius2);
604  jetCont2->SetRhoName(nrho2);
605  jetCont2->SetLeadingHadronType(biasType);
606  jetCont2->SetPtBiasJetTrack(jetBias);
607  jetCont2->SetPtBiasJetClus(jetBias);
608  jetCont2->SetJetPtCut(jetptcut);
609  jetCont2->SetPercAreaCut(jetareacut);
610  jetCont2->SetIsParticleLevel(kTRUE);
611  jetCont2->ConnectParticleContainer(trackCont2);
612  jetCont2->ConnectClusterContainer(clusCont2);
613  jetCont2->SetMaxTrackPt(1000); // disable default 100 GeV/c track cut for particle level jets
614 
615  jetTask->SetMatching(matching, maxDistance1, maxDistance2);
616  jetTask->SetVzRange(-10,10);
617  jetTask->SetIsPythia(kTRUE);
618  jetTask->SetPtHardBin(ptHardBin);
619  jetTask->SetCentRange(minCent,maxCent);
620 
621  if (biggerMatrix) jetTask->SetHistoBins(1000,0,500);
622 
623  //-------------------------------------------------------
624  // Final settings, pass to manager and set the containers
625  //-------------------------------------------------------
626 
627  mgr->AddTask(jetTask);
628 
629  // Create containers for input/output
630  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
631  TString contname(name);
632  contname += "_histos";
633  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
634  TList::Class(),AliAnalysisManager::kOutputContainer,
635  Form("%s", AliAnalysisManager::GetCommonFileName()));
636  mgr->ConnectInput (jetTask, 0, cinput1 );
637  mgr->ConnectOutput (jetTask, 1, coutput1 );
638 
639  return jetTask;
640 }
641 
642 void AliAnalysisTaskSoftDropResponse::Decluster(const fastjet::PseudoJet& jet) {
643 
644  fastjet::PseudoJet jet1;
645  fastjet::PseudoJet jet2;
646 
647  if ( jet.has_parents(jet1, jet2) ) {
648 
649  ++fNsdsteps;
650 
651  Float_t pt1 = jet1.pt();
652  Float_t pt2 = jet2.pt();
653 
654  Float_t dr = TMath::Sqrt( jet1.plain_distance(jet2) );
655 
656  Float_t z;
657  if (pt1 < pt2) z = pt1/(pt1+pt2);
658  else z = pt2/(pt1+pt2);
659 
660  if (z > 0.1) {
661  fZg_values.push_back(z);
662  fRg_values.push_back(dr);
663  fSDsteps_values.push_back(fNsdsteps);
664  }
665 
666  if (pt1 > pt2) Decluster(jet1);
667  else Decluster(jet2);
668 
669  }
670 
671 }
672 
673 //________________________________________________________________________
674 fastjet::ClusterSequence* AliAnalysisTaskSoftDropResponse::Recluster(const AliEmcalJet* jet) {
675 
676  std::vector<fastjet::PseudoJet> particles;
677  UShort_t ntracks = jet->GetNumberOfTracks();
678  for (int j = 0; j < ntracks; j++) {
679  particles.push_back( fastjet::PseudoJet( jet->Track(j)->Px(), jet->Track(j)->Py(), jet->Track(j)->Pz(), jet->Track(j)->E() ) );
680  }
681  fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, 0.4, fastjet::E_scheme);
682  return ( new fastjet::ClusterSequence(particles, jet_def) );
683 }
684 
686  fNsdsteps = 0;
687  fZg_values.clear();
688  fRg_values.clear();
689  fSDsteps_values.clear();
690  fastjet::ClusterSequence* cs = Recluster(jet);
691  if (cs) {
692  std::vector<fastjet::PseudoJet> jetrecl = sorted_by_pt( cs->inclusive_jets() );
693  if (jetrecl.size() > 0) Decluster( jetrecl[0] );
694  }
695  delete cs;
696 }
697 
698 //________________________________________________________________________
700 {
701  // Create user objects.
702 
704 
705  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
706  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
707 
708  if (!jets1 || !jets2) return;
709 
710  if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
711  else fIsJet1Rho = kTRUE;
712  if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
713  else fIsJet2Rho = kTRUE;
714 
715  fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250);
716  fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason");
717  fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
718  fHistRejectionReason1->GetZaxis()->SetTitle("counts");
721 
722  fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250);
723  fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason");
724  fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
725  fHistRejectionReason2->GetZaxis()->SetTitle("counts");
728 
729  if (fHistoType==0)
730  AllocateTH2();
731  else
733 
734  // Initialize
736  if (embeddingHelper) {
737  bool res = fEmbeddingQA.Initialize();
738  if (res) {
740  }
741  }
742 
743  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
744 }
745 
746 
747 //________________________________________________________________________
749 {
750  // Destructor
751 }
void SetCentRange(Double_t min, Double_t max)
Double_t Area() const
Definition: AliEmcalJet.h:130
TH2 * fHistRejectionReason1
whether the jet2 collection has to be average subtracted
void SetParticlePtCut(Double_t cut)
Bool_t FillHistograms()
Function filling histograms.
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
Double_t MCPt() const
Definition: AliEmcalJet.h:153
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:328
Double_t Eta() const
Definition: AliEmcalJet.h:121
void SetLeadingHadronType(Int_t t)
void FillJetHisto(AliEmcalJet *jet, Int_t Set)
Double_t Phi() const
Definition: AliEmcalJet.h:117
void SetPtBiasJetTrack(Float_t b)
AliEmcalEmbeddingQA fEmbeddingQA
! Embedding QA hists (will only be added if embedding)
AliEmcalJet * MatchedJet() const
Definition: AliEmcalJet.h:331
void GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
std::vector< Float_t > fRg_values
! n_SD ordering
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
AliVParticle * Track(Int_t idx) const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
THnSparse * fHistMatching
jet2 THnSparse
void SetPercAreaCut(Float_t p)
void SetVzRange(Double_t min, Double_t max)
void SetPtBiasJetClus(Float_t b)
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
Bool_t fIsEmbedded
trigger, embedded signal
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Bool_t fIsJet2Rho
whether the jet1 collection has to be average subtracted
bool AddQAPlotsToList(TList *list)
void SetRhoName(const char *n)
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
std::vector< Float_t > fZg_values
! n_SD ordering
unsigned int UInt_t
Definition: External.C:33
THnSparse * fHistJets1
Rejection reason vs. jet pt.
float Float_t
Definition: External.C:68
UShort_t GetMatchingType() const
Definition: AliEmcalJet.h:332
void Decluster(const fastjet::PseudoJet &jet)
void SetPtHardBin(Int_t b)
Implementation of task to embed external events.
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
void SetMatching(MatchingType t, Double_t p1=1, Double_t p2=1)
static AliAnalysisTaskSoftDropResponse * AddTaskSoftDropResponse(const char *ntracks1="Tracks", const char *nclusters1="CaloClusters", const char *njets1="Jets", const char *nrho1="Rho", const Double_t jetradius1=0.4, const char *ntracks2="MCParticles", const char *nclusters2="", const char *njets2="MCJets", const char *nrho2="", const Double_t jetradius2=0.4, const Double_t jetptcut=1, const Double_t jetareacut=0.557, const Double_t jetBias=5, const Int_t biasType=0, const AliAnalysisTaskSoftDropResponse::MatchingType matching=AliAnalysisTaskSoftDropResponse::kGeometrical, const Double_t maxDistance1=0.25, const Double_t maxDistance2=0.25, const char *cutType="TPC", const Int_t ptHardBin=-999, const Double_t minCent=-999, const Double_t maxCent=-999, const char *taskname="AliAnalysisTaskSoftDropResponse", const Bool_t biggerMatrix=kFALSE, AliAnalysisTaskSoftDropResponse *address=0, const Double_t maxTrackPt=100)
Double_t GetJetPtCut() const
TH2 * fHistRejectionReason2
Rejection reason vs. jet pt.
TObjArray fJetCollArray
jet collection array
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalJet * GetNextJet()
AliEmcalList * fOutput
!output list
void SetHistoBins(Int_t nbins, Double_t min, Double_t max)
void GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
THnSparse * fHistJets2
jet1 THnSparse
void SetMakeGeneralHistograms(Bool_t g)
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
unsigned short UShort_t
Definition: External.C:28
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
void ConnectClusterContainer(AliClusterContainer *c)
void SetMaxTrackPt(Float_t b)
Int_t fZ2gAxis
add Z2g axis in matching THnSparse (default=0)
Container structure for EMCAL clusters.
void GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
void FillMatchingHistos(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t d, Double_t CE1, Double_t CE2)
fastjet::ClusterSequence * Recluster(const AliEmcalJet *jet)
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()