AliPhysics  b555aef (b555aef)
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"
25 
27 
28 //________________________________________________________________________
31  fZ2gAxis(0)
32 {
33  // Default constructor.
34 
35  SetMakeGeneralHistograms(kTRUE);
36 }
37 
38 //________________________________________________________________________
40  AliJetResponseMaker(name),
41  fZ2gAxis(0)
42 {
43  // Standard constructor.
44 
46 }
47 
48 
49 //________________________________________________________________________
51 {
52  // Allocate THnSparse histograms.
53 
54  TString title[25]= {""};
55  Int_t nbins[25] = {0};
56  Double_t min[25] = {0.};
57  Double_t max[25] = {0.};
58  Int_t dim = 0;
59 
60  title[dim] = "#phi";
61  nbins[dim] = fNbins/4;
62  min[dim] = 0;
63  max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
64  dim++;
65 
66  title[dim] = "#eta";
67  nbins[dim] = fNbins/4;
68  min[dim] = -1;
69  max[dim] = 1;
70  dim++;
71 
72  title[dim] = "p_{T}";
73  nbins[dim] = fNbins;
74  min[dim] = 0;
75  max[dim] = 250;
76  dim++;
77 
78  title[dim] = "A_{jet}";
79  nbins[dim] = fNbins/4;
80  min[dim] = 0;
81  max[dim] = 1.5;
82  dim++;
83 
84  Int_t dim1 = dim, dim2 = dim;
85 
86  if (fIsJet1Rho) {
87  title[dim1] = "p_{T}^{corr}";
88  nbins[dim1] = fNbins*2;
89  min[dim1] = -250;
90  max[dim1] = 250;
91  dim1++;
92  }
93 
94  if (fIsEmbedded) {
95  title[dim1] = "p_{T}^{MC}";
96  nbins[dim1] = fNbins;
97  min[dim1] = 0;
98  max[dim1] = 250;
99  dim1++;
100  }
101 
102  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
103  for (Int_t i = 0; i < dim1; i++) fHistJets1->GetAxis(i)->SetTitle(title[i]);
104  fOutput->Add(fHistJets1);
105 
106  if (fIsJet2Rho) {
107  title[dim2] = "p_{T}^{corr}";
108  nbins[dim2] = fNbins*2;
109  min[dim2] = -250;
110  max[dim2] = 250;
111  dim2++;
112  }
113 
114  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
115  for (Int_t i = 0; i < dim2; i++) fHistJets2->GetAxis(i)->SetTitle(title[i]);
116  fOutput->Add(fHistJets2);
117 
118  // Matching
119 
120  dim = 0;
121 
122  title[dim] = "p_{T,1}";
123  nbins[dim] = fNbins;
124  min[dim] = 0;
125  max[dim] = 250;
126  dim++;
127 
128  title[dim] = "p_{T,2}";
129  nbins[dim] = fNbins;
130  min[dim] = 0;
131  max[dim] = 250;
132  dim++;
133 
134  title[dim] = "A_{jet,1}";
135  nbins[dim] = fNbins/4;
136  min[dim] = 0;
137  max[dim] = 1.5;
138  dim++;
139 
140  title[dim] = "A_{jet,2}";
141  nbins[dim] = fNbins/4;
142  min[dim] = 0;
143  max[dim] = 1.5;
144  dim++;
145 
146  title[dim] = "distance";
147  nbins[dim] = fNbins/2;
148  min[dim] = 0;
149  max[dim] = 1.2;
150  dim++;
151 
152  title[dim] = "CE1";
153  nbins[dim] = fNbins/2;
154  min[dim] = 0;
155  max[dim] = 1.2;
156  dim++;
157 
158  title[dim] = "CE2";
159  nbins[dim] = fNbins/2;
160  min[dim] = 0;
161  max[dim] = 1.2;
162  dim++;
163 
164  if (fDeltaPtAxis) {
165  title[dim] = "#deltaA_{jet}";
166  nbins[dim] = fNbins/2;
167  min[dim] = -1;
168  max[dim] = 1;
169  dim++;
170 
171  title[dim] = "#deltap_{T}";
172  nbins[dim] = fNbins*2;
173  min[dim] = -250;
174  max[dim] = 250;
175  dim++;
176  }
177  if (fIsJet1Rho) {
178  title[dim] = "p_{T,1}^{corr}";
179  nbins[dim] = fNbins*2;
180  min[dim] = -250;
181  max[dim] = 250;
182  dim++;
183  }
184  if (fIsJet2Rho) {
185  title[dim] = "p_{T,2}^{corr}";
186  nbins[dim] = fNbins*2;
187  min[dim] = -250;
188  max[dim] = 250;
189  dim++;
190  }
191  if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
192  title[dim] = "#deltap_{T}^{corr}";
193  nbins[dim] = fNbins*2;
194  min[dim] = -250;
195  max[dim] = 250;
196  dim++;
197  }
198  if (fDeltaEtaDeltaPhiAxis) {
199  title[dim] = "#delta#eta";
200  nbins[dim] = fNbins/2;
201  min[dim] = -1;
202  max[dim] = 1;
203  dim++;
204 
205  title[dim] = "#delta#phi";
206  nbins[dim] = fNbins/2;
207  min[dim] = -TMath::Pi()/2;
208  max[dim] = TMath::Pi()*3/2;
209  dim++;
210  }
211  if (fIsEmbedded) {
212  title[dim] = "p_{T,1}^{MC}";
213  nbins[dim] = fNbins;
214  min[dim] = 0;
215  max[dim] = 250;
216  dim++;
217 
218  if (fDeltaPtAxis) {
219  title[dim] = "#deltap_{T}^{MC}";
220  nbins[dim] = fNbins*2;
221  min[dim] = -250;
222  max[dim] = 250;
223  dim++;
224  }
225  }
226 
227  if (fZgAxis) {
228  title[dim] = "Z_{g,1}";
229  nbins[dim] = 20;
230  min[dim] = 0.0;
231  max[dim] = 1.0;
232  dim++;
233  title[dim] = "Z_{g,2}";
234  nbins[dim] = 20;
235  min[dim] = 0.0;
236  max[dim] = 1.0;
237  dim++;
238  }
239 
240  if (fZ2gAxis) {
241  title[dim] = "Z2_{g,1}";
242  nbins[dim] = 20;
243  min[dim] = 0.0;
244  max[dim] = 1.0;
245  dim++;
246  title[dim] = "Z2_{g,2}";
247  nbins[dim] = 20;
248  min[dim] = 0.0;
249  max[dim] = 1.0;
250  dim++;
251  }
252 
253  if (fdRAxis) {
254  title[dim] = "dR_{1}";
255  nbins[dim] = 40;
256  min[dim] = 0.0;
257  max[dim] = 0.5;
258  dim++;
259  title[dim] = "dR_{2}";
260  nbins[dim] = 40;
261  min[dim] = 0.0;
262  max[dim] = 0.5;
263  dim++;
264  }
265 
266  if (fPtgAxis) {
267  title[dim] = "p_{T,g,1}";
268  nbins[dim] = 16;
269  min[dim] = 0.0;
270  max[dim] = 160.0;
271  dim++;
272  title[dim] = "p_{T,g,2}";
273  nbins[dim] = 16;
274  min[dim] = 0.0;
275  max[dim] = 160.0;
276  dim++;
277  }
278 
279  if (fDBCAxis) {
280  title[dim] = "DBC_{1}";
281  nbins[dim] = 20;
282  min[dim] = 0.0;
283  max[dim] = 20.0;
284  dim++;
285  title[dim] = "DBC_{2}";
286  nbins[dim] = 20;
287  min[dim] = 0.0;
288  max[dim] = 20.0;
289  dim++;
290  }
291 
292  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
293 
294  for (Int_t i = 0; i < dim; i++)
295  fHistMatching->GetAxis(i)->SetTitle(title[i]);
296 
297  fOutput->Add(fHistMatching);
298 }
299 
300 //________________________________________________________________________
302 {
303 
304  zg = -1.0;
305 
306  std::vector<fastjet::PseudoJet> particles;
307  UShort_t ntracks = jet->GetNumberOfTracks();
308  for (int j = 0; j < ntracks; j++) {
309  particles.push_back( fastjet::PseudoJet( jet->Track(j)->Px(), jet->Track(j)->Py(), jet->Track(j)->Pz(), jet->Track(j)->E() ) );
310  }
311  fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, 0.4, fastjet::E_scheme);
312  fastjet::ClusterSequence cs(particles, jet_def);
313  std::vector<fastjet::PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
314 
315  if (jets.size() > 0) {
316  zg = AliAnalysisTaskSoftDrop::SoftDropDeclustering(jets[0], zcut, beta);
317  }
318 
319 }
320 
321 //________________________________________________________________________
323 {
324  AliJetContainer* jets1 = GetJetContainer(0);
325  AliJetContainer* jets2 = GetJetContainer(1);
326 
327  Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
328  Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
329 
330  Float_t z2g1 = -1.0;
331  Float_t z2g2 = -1.0;
332  CalculateZg(jet1, 0.5, 1.5, z2g1);
333  CalculateZg(jet2, 0.5, 1.5, z2g2);
334 
335  if (fHistoType==1) {
336  Double_t contents[20]={0};
337 
338  for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
339  TString title(fHistMatching->GetAxis(i)->GetTitle());
340  if (title=="p_{T,1}")
341  contents[i] = jet1->Pt();
342  else if (title=="p_{T,2}")
343  contents[i] = jet2->Pt();
344  else if (title=="A_{jet,1}")
345  contents[i] = jet1->Area();
346  else if (title=="A_{jet,2}")
347  contents[i] = jet2->Area();
348  else if (title=="distance")
349  contents[i] = d;
350  else if (title=="CE1")
351  contents[i] = CE1;
352  else if (title=="CE2")
353  contents[i] = CE2;
354  else if (title=="#deltaA_{jet}")
355  contents[i] = jet1->Area()-jet2->Area();
356  else if (title=="#deltap_{T}")
357  contents[i] = jet1->Pt()-jet2->Pt();
358  else if (title=="#delta#eta")
359  contents[i] = jet1->Eta()-jet2->Eta();
360  else if (title=="#delta#phi")
361  contents[i] = jet1->Phi()-jet2->Phi();
362  else if (title=="p_{T,1}^{corr}")
363  contents[i] = corrpt1;
364  else if (title=="p_{T,2}^{corr}")
365  contents[i] = corrpt2;
366  else if (title=="#deltap_{T}^{corr}")
367  contents[i] = corrpt1-corrpt2;
368  else if (title=="p_{T,1}^{MC}")
369  contents[i] = jet1->MCPt();
370  else if (title=="#deltap_{T}^{MC}")
371  contents[i] = jet1->MCPt()-jet2->Pt();
372  else if (title=="Z_{g,1}")
373  contents[i] = jet1->GetShapeProperties()->GetSoftDropZg();
374  else if (title=="Z_{g,2}")
375  contents[i] = jet2->GetShapeProperties()->GetSoftDropZg();
376  else if (title=="Z2_{g,1}")
377  contents[i] = z2g1;
378  else if (title=="Z2_{g,2}")
379  contents[i] = z2g2;
380  else if (title=="dR_{1}")
381  contents[i] = jet1->GetShapeProperties()->GetSoftDropdR();
382  else if (title=="dR_{2}")
383  contents[i] = jet2->GetShapeProperties()->GetSoftDropdR();
384  else if (title=="p_{T,g,1}")
385  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropPtfrac() )*( jet1->Pt() );
386  else if (title=="p_{T,g,2}")
387  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropPtfrac() )*( jet2->Pt() );
388  else if (title=="DBC_{1}")
389  contents[i] = ( jet1->GetShapeProperties()->GetSoftDropDropCount() );
390  else if (title=="DBC_{2}")
391  contents[i] = ( jet2->GetShapeProperties()->GetSoftDropDropCount() );
392  else
393  AliWarning(Form("Unable to fill dimension %s!",title.Data()));
394  }
395 
396  fHistMatching->Fill(contents);
397  }
398 
399 }
400 
405  const char *ntracks1,
406  const char *nclusters1,
407  const char *njets1,
408  const char *nrho1,
409  const Double_t jetradius1,
410  const char *ntracks2,
411  const char *nclusters2,
412  const char *njets2,
413  const char *nrho2,
414  const Double_t jetradius2,
415  const Double_t jetptcut,
416  const Double_t jetareacut,
417  const Double_t jetBias,
418  const Int_t biasType,
420  const Double_t maxDistance1,
421  const Double_t maxDistance2,
422  const char *cutType,
423  const Int_t ptHardBin,
424  const Double_t minCent,
425  const Double_t maxCent,
426  const char *taskname,
427  const Bool_t biggerMatrix,
429  const Double_t maxTrackPt)
430 {
431  // Get the pointer to the existing analysis manager via the static access method.
432  //==============================================================================
433  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
434  if (!mgr)
435  {
436  ::Error("AddTaskJetResponseMaker", "No analysis manager to connect to.");
437  return NULL;
438  }
439 
440  // Check the analysis type using the event handlers connected to the analysis manager.
441  //==============================================================================
442  if (!mgr->GetInputEventHandler())
443  {
444  ::Error("AddTaskJetResponseMaker", "This task requires an input event handler");
445  return NULL;
446  }
447 
448  //-------------------------------------------------------
449  // Init the task and do settings
450  //-------------------------------------------------------
451 
452  TString name(Form("%s_%s_%s_Bias%d_BiasType%d_%s",taskname,njets1,njets2,(Int_t)floor(jetBias),biasType,cutType));
453 
454  if (minCent != -999 && maxCent != -999)
455  name += Form("_Cent%d_%d", (Int_t)floor(minCent), (Int_t)floor(maxCent));
456 
457  if (ptHardBin != -999)
458  name += Form("_PtHard%d", ptHardBin);
459 
460  AliAnalysisTaskSoftDropResponse* jetTask = address;
461  if (jetTask) {
462  new (jetTask) AliAnalysisTaskSoftDropResponse(name);
463  }
464  else {
465  jetTask = new AliAnalysisTaskSoftDropResponse(name);
466  }
467 
468  AliParticleContainer *trackCont1 = jetTask->AddParticleContainer(ntracks1);
469  AliClusterContainer *clusCont1 = jetTask->AddClusterContainer(nclusters1);
470  AliJetContainer *jetCont1 = jetTask->AddJetContainer(njets1, cutType, jetradius1);
471  jetCont1->SetRhoName(nrho1);
472  jetCont1->SetLeadingHadronType(biasType);
473  jetCont1->SetPtBiasJetTrack(jetBias);
474  jetCont1->SetPtBiasJetClus(jetBias);
475  jetCont1->SetJetPtCut(jetptcut);
476  jetCont1->SetPercAreaCut(jetareacut);
477  jetCont1->SetIsParticleLevel(kFALSE);
478  jetCont1->ConnectParticleContainer(trackCont1);
479  jetCont1->ConnectClusterContainer(clusCont1);
480  jetCont1->SetMaxTrackPt(maxTrackPt);
481 
482  AliParticleContainer *trackCont2 = jetTask->AddParticleContainer(ntracks2);
483  trackCont2->SetParticlePtCut(0);
484  AliClusterContainer *clusCont2 = jetTask->AddClusterContainer(nclusters2);
485  AliJetContainer *jetCont2 = jetTask->AddJetContainer(njets2, cutType, jetradius2);
486  jetCont2->SetRhoName(nrho2);
487  jetCont2->SetLeadingHadronType(biasType);
488  jetCont2->SetPtBiasJetTrack(jetBias);
489  jetCont2->SetPtBiasJetClus(jetBias);
490  jetCont2->SetJetPtCut(jetptcut);
491  jetCont2->SetPercAreaCut(jetareacut);
492  jetCont2->SetIsParticleLevel(kTRUE);
493  jetCont2->ConnectParticleContainer(trackCont2);
494  jetCont2->ConnectClusterContainer(clusCont2);
495  jetCont2->SetMaxTrackPt(1000); // disable default 100 GeV/c track cut for particle level jets
496 
497  jetTask->SetMatching(matching, maxDistance1, maxDistance2);
498  jetTask->SetVzRange(-10,10);
499  jetTask->SetIsPythia(kTRUE);
500  jetTask->SetPtHardBin(ptHardBin);
501  jetTask->SetCentRange(minCent,maxCent);
502 
503  if (biggerMatrix) jetTask->SetHistoBins(1000,0,500);
504 
505  //-------------------------------------------------------
506  // Final settings, pass to manager and set the containers
507  //-------------------------------------------------------
508 
509  mgr->AddTask(jetTask);
510 
511  // Create containers for input/output
512  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
513  TString contname(name);
514  contname += "_histos";
515  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname.Data(),
516  TList::Class(),AliAnalysisManager::kOutputContainer,
517  Form("%s", AliAnalysisManager::GetCommonFileName()));
518  mgr->ConnectInput (jetTask, 0, cinput1 );
519  mgr->ConnectOutput (jetTask, 1, coutput1 );
520 
521  return jetTask;
522 }
523 
524 
525 //________________________________________________________________________
527 {
528  // Destructor
529 }
void SetCentRange(Double_t min, Double_t max)
Double_t Area() const
Definition: AliEmcalJet.h:130
void SetParticlePtCut(Double_t cut)
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:27
Double_t MCPt() const
Definition: AliEmcalJet.h:153
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
void SetLeadingHadronType(Int_t t)
Double_t Phi() const
Definition: AliEmcalJet.h:117
void SetPtBiasJetTrack(Float_t b)
AliVParticle * Track(Int_t idx) 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
void SetRhoName(const char *n)
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
THnSparse * fHistJets1
Rejection reason vs. jet pt.
float Float_t
Definition: External.C:68
void CalculateZg(AliEmcalJet *jet, const Float_t zcut, const Float_t beta, Float_t zg)
void SetPtHardBin(Int_t b)
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)
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalList * fOutput
!output list
void SetHistoBins(Int_t nbins, Double_t min, Double_t max)
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
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)
static Float_t SoftDropDeclustering(fastjet::PseudoJet jet, const Float_t zcut, const Float_t beta)
Container structure for EMCAL clusters.
void FillMatchingHistos(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t d, Double_t CE1, Double_t CE2)
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins