AliPhysics  master (3d17d9d)
AliAnalysisTaskSDKLResponse.cxx
Go to the documentation of this file.
1 #include <TClonesArray.h>
2 #include <TMath.h>
3 #include <TH1F.h>
4 #include <TH2F.h>
5 #include <TH3F.h>
6 #include <THnSparse.h>
7 #include <TNtuple.h>
8 #include <TList.h>
9 #include <TLorentzVector.h>
10 #include <TRandom.h>
11 
12 #include "AliVCluster.h"
13 //#include "AliAODCaloCluster.h"
14 //#include "AliESDCaloCluster.h"
15 #include "AliVTrack.h"
16 #include "AliAODMCParticle.h"
17 #include "AliEmcalJet.h"
18 #include "AliRhoParameter.h"
19 #include "AliLog.h"
20 #include "AliJetContainer.h"
21 #include "AliParticleContainer.h"
22 //#include "AliClusterContainer.h"
23 //#include "AliPicoTrack.h"
24 
25 #include "AliAnalysisManager.h"
26 #include "AliAnalysisTaskSDKL.h"
28 
29 //ClassImp(AliAnalysisTaskSDKLResponse)
30 
31 //________________________________________________________________________
33  AliAnalysisTaskSDKL("AliAnalysisTaskSDKLResponse"),
34  fhResponse(),
35  fhPtDeltaPt(),
36  fhPtDeltaZg(),
37  fhPtDeltaRg(),
38  fhPtDeltaMg(),
39  fhResponseBackSub(),
40  fhPtDeltaPtBackSub(0),
41  fhPtDeltaZgBackSub(),
42  fhPtDeltaRgBackSub(),
43  fhPtDeltaMgBackSub(),
44  fhResponseDet(),
45  fhPtDeltaPtDet(0),
46  fhPtDeltaZgDet(),
47  fhPtDeltaRgDet(),
48  fhPtDeltaMgDet(),
49 // fhRho(0),
50 // fhRhoSparse(0),
51  fhPtDeltaPtAreaBackSub(0),
52  fhPtDeltaPtAreaBackSubSparse(0),
53  fTreeDL(0),
54  fTreeDLUEBS(0),
55  fTreePL(0),
56  fJetsCont1(0),
57  fJetsCont2(0),
58  fTracksCont1(0),
59  fTracksCont2(0),
60  fFractionEventsDumpedToTree(0),
61  fRandom(0)
62 {
63  // Default constructor.
64  // SetMakeGeneralHistograms(kTRUE);
65 }
66 
67 //________________________________________________________________________
68 AliAnalysisTaskSDKLResponse::AliAnalysisTaskSDKLResponse(const char *name, Int_t const backgroption, Double_t const fractioneventsfortree) :
69  AliAnalysisTaskSDKL(name, backgroption),
70  fhResponse(),
71  fhPtDeltaPt(0),
72  fhPtDeltaZg(),
73  fhPtDeltaRg(),
74  fhPtDeltaMg(),
80  fhResponseDet(),
81  fhPtDeltaPtDet(0),
85 // fhRho(0),
86 // fhRhoSparse(0),
89  fTreeDL(0),
90  fTreeDLUEBS(0),
91  fTreePL(0),
92  fJetsCont1(0),
93  fJetsCont2(0),
94  fTracksCont1(0),
95  fTracksCont2(0),
96  fFractionEventsDumpedToTree(fractioneventsfortree),
97  fRandom(0)
98 {
99  // Standard constructor.
100  // SetMakeGeneralHistograms(kTRUE);
101  DefineOutput(2, TTree::Class());
102  DefineOutput(3, TTree::Class());
103  DefineOutput(4, TTree::Class());
104 }
105 
106 //________________________________________________________________________
108 {
109  if (fTreeDL) delete fTreeDL;
110  if (fTreeDLUEBS) delete fTreeDLUEBS;
111  if (fTreePL) delete fTreePL;
112  if (fRandom) delete fRandom;
113  // Destructor.
114 }
115 
117  const char *ntracks,
118  const char *njets1,
119  const char *njets2,
120  const char *nrho,
122  Double_t jetradius,
123  Double_t jetptcut,
124  Double_t jetareacut,
125  const char *type,
126  Int_t backgroption,
127  Int_t leadhadtype,
128  Double_t fractioneventsfortree,
129  const char *taskname
130 )
131 {
132  // Get the pointer to the existing analysis manager via the static access method.
133  //==============================================================================
134  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
135  if (!mgr)
136  {
137  ::Error("AddTaskEmcalJetSample", "No analysis manager to connect to.");
138  return NULL;
139  }
140 
141  // Check the analysis type using the event handlers connected to the analysis manager.
142  //==============================================================================
143  if (!mgr->GetInputEventHandler())
144  {
145  ::Error("AddTaskEmcalJetSample", "This task requires an input event handler");
146  return NULL;
147  }
148 
149  //-------------------------------------------------------
150  // Init the task and do settings
151  //-------------------------------------------------------
152 
153  TString trackName(ntracks);
154 
155  TString name(taskname);
156  if (strcmp(njets1,"")) {
157  name += "_";
158  name += njets1;
159  }
160  if (strcmp(njets2,"")) {
161  name += "_";
162  name += njets2;
163  }
164  if (strcmp(nrho,"")) {
165  name += "_";
166  name += nrho;
167  }
168  if (strcmp(type,"")) {
169  name += "_";
170  name += type;
171  }
172 
173  Printf("name: %s",name.Data());
174 
175  AliAnalysisTaskSDKLResponse* jetTask = new AliAnalysisTaskSDKLResponse(name, backgroption, fractioneventsfortree);
176  //jetTask->SetCentRange(0.,100.);
177  //jetTask->SetNCentBins(nCentBins);
178  //AliMultSelection *multSelection =static_cast<AliMultSelection*>(fAOD->FindListObject("MultSelection"));
179  //if(multSelection) centrality = multSelection->GetMultiplicityPercentile("V0M");
180 
181 // AliParticleContainer* partCont = 0;
182 // if (trackName == "mcparticles") {
183 // partCont = jetTask->AddMCParticleContainer(ntracks);
184 // }
185 // else if (trackName == "usedefault" || trackName == "tracks" || trackName == "Tracks") {
186 // partCont = jetTask->AddTrackContainer(ntracks);
187 // }
188 // else if (!trackName.IsNull()) {
189 // partCont = new AliParticleContainer(trackName);
190 // }
191 
192 // AliParticleContainer* partCont_pl = jetTask->AddParticleContainer("mcparticles");
193  AliParticleContainer* partCont = jetTask->AddTrackContainer("tracks");
194  AliParticleContainer* partCont_emb = jetTask->AddTrackContainer("tracks");
195 
196  TString strType(type);
197  AliJetContainer *jetCont1 = jetTask->AddJetContainer(njets1,strType,jetradius);
198  if(jetCont1) {
199  jetCont1->SetRhoName(nrho);
200 // jetCont1->ConnectParticleContainer(partCont);
201  //jetCont->SetZLeadingCut(0.98,0.98);
202  jetCont1->SetPercAreaCut(jetareacut);
203  jetCont1->SetJetPtCut(jetptcut);
204  jetCont1->SetLeadingHadronType(leadhadtype);
205  }
206 
207  AliJetContainer *jetCont2 = jetTask->AddJetContainer(njets2,strType,jetradius);
208  if(jetCont2) {
209  jetCont2->SetRhoName(nrho);
210  //jetCont2->ConnectParticleContainer(partCont);
211  jetCont2->SetPercAreaCut(jetareacut);
212  jetCont2->SetJetPtCut(jetptcut);
213  jetCont2->SetLeadingHadronType(leadhadtype);
214  }
215 
216  //-------------------------------------------------------
217  // Final settings, pass to manager and set the containers
218  //-------------------------------------------------------
219 
220  mgr->AddTask(jetTask);
221 
222  // Create containers for input/output
223  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer() ;
224  TString contname1(name);
225  contname1 += "_histos";
226  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contname1.Data(),
227  TList::Class(),AliAnalysisManager::kOutputContainer,
228  Form("%s", AliAnalysisManager::GetCommonFileName()));
229 
230  TString contname2(name);
231  contname2 += "_treeDL";
232  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(contname2.Data(),
233  TTree::Class(),AliAnalysisManager::kOutputContainer,
234  Form("%s", AliAnalysisManager::GetCommonFileName()));
235 
236  TString contname3(name);
237  contname3 += "_treeDLUEBS";
238  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(contname3.Data(),
239  TTree::Class(),AliAnalysisManager::kOutputContainer,
240  Form("%s", AliAnalysisManager::GetCommonFileName()));
241 
242  TString contname4(name);
243  contname4 += "_treePL";
244  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(contname4.Data(),
245  TTree::Class(),AliAnalysisManager::kOutputContainer,
246  Form("%s", AliAnalysisManager::GetCommonFileName()));
247 
248  mgr->ConnectInput (jetTask, 0, cinput1 );
249  mgr->ConnectOutput (jetTask, 1, coutput1 );
250  mgr->ConnectOutput (jetTask, 2, coutput2 );
251  mgr->ConnectOutput (jetTask, 3, coutput3 );
252  mgr->ConnectOutput (jetTask, 4, coutput4 );
253 
254  return jetTask;
255 }
256 
257 //________________________________________________________________________
259  // Create user output.
260 
262  OpenFile(1);
263 
264  const Int_t nbins = 11;
265 
266  //dist eshare iscl pt1 zg1 rg1 mg1/pt1 nsd1 pt2 zg2 rg2 mg2/pt2
267  // 0 1 2 3 4 5 6 7 8 9 10
268  Int_t bins[nbins] = { 10, 10, 2, 200, 20, 20, 20, 200, 20, 20, 20};
269  Double_t xmin[nbins] = {0.0, 0.0, 0, 0.0, 0.0, 0.0, 0.0, 0., 0.0, 0.0, 0.0};
270  Double_t xmax[nbins] = {0.8, 1.0, 2, 200., 0.5, 0.4, 0.2, 200., 0.5, 0.4, 0.2};
271 
272  fhResponse[0] = new THnSparseD("hResponse_1", "hResponse_1", nbins, bins, xmin, xmax); //by splits
273  fhResponse[1] = new THnSparseD("hResponse_2", "hResponse_2", nbins, bins, xmin, xmax);
274 
275  fhResponseBackSub[0] = new THnSparseD("hResponseBackSub_1", "hResponseBackSub_1", nbins, bins, xmin, xmax);
276  fhResponseBackSub[1] = new THnSparseD("hResponseBackSub_2", "hResponseBackSub_2", nbins, bins, xmin, xmax);
277 
278  fhResponseDet[0] = new THnSparseD("hResponseDet_1", "hResponseDet_1", nbins, bins, xmin, xmax);
279  fhResponseDet[1] = new THnSparseD("hResponseDet_2", "hResponseDet_2", nbins, bins, xmin, xmax);
280 
281  xmax[5] = 0.3; //rg1
282  xmax[9] = 0.3; //rg2
283  xmax[6] = 0.15; //mg1/pt
284  xmax[10] = 0.15; //mg2/pt
285  fhResponse[2] = new THnSparseD("hResponse_3", "hResponse_3", nbins, bins, xmin, xmax);
286  fhResponseBackSub[2] = new THnSparseD("hResponseBackSub_3", "hResponseBackSub_3", nbins, bins, xmin, xmax);
287  fhResponseDet[2] = new THnSparseD("hResponseDet_3", "hResponseDet_3", nbins, bins, xmin, xmax);
288 
289  xmax[6] = 0.1; //mg1/pt
290  xmax[10] = 0.1; //mg2/pt
291  fhResponse[3] = new THnSparseD("hResponse_4", "hResponse_4", nbins, bins, xmin, xmax);
292  fhResponseBackSub[3] = new THnSparseD("hResponseBackSub_4", "hResponseBackSub_4", nbins, bins, xmin, xmax);
293  fhResponseDet[3] = new THnSparseD("hResponseDet_4", "hResponseDet_4", nbins, bins, xmin, xmax);
294 
295  for (int s = 0; s < 4; s++) {
296  fOutput->Add( fhResponse[s] );
297  fOutput->Add( fhResponseBackSub[s] );
298  fOutput->Add( fhResponseDet[s] );
299  }
300 
301  fhPtDeltaPt = new TH2F("fhPtDeltaPt","fhPtDeltaPt",15,0.,150.,200,-20.,20.);
302  fOutput->Add(fhPtDeltaPt);
303 
304  fhPtDeltaPtBackSub = new TH2F("fhPtDeltaPtBackSub","fhPtDeltaPtBackSub",15,0.,150.,200,-20.,20.);
306 
307  fhPtDeltaPtDet = new TH2F("fhPtDeltaPtDet","fhPtDeltaPtDet",15,0.,150.,200,-20.,20.);
308  fOutput->Add(fhPtDeltaPtDet);
309 
310  auto r_l = -0.05;
311  auto r_r = 0.05;
312  for (int s = 0; s < 4; s++) {
313 
314  TString dzg_name = "fhPtDeltaZg";
315  dzg_name += (s+1);
316  fhPtDeltaZg[s] = new TH2F(dzg_name,dzg_name,15,0.,150.,200,-0.2,0.2);
317 
318  TString drg_name = "fhPtDeltaRg";
319  drg_name += (s+1);
320  fhPtDeltaRg[s] = new TH2F(drg_name,drg_name,15,0.,150.,200,r_l,r_r);
321 
322  TString dmg_name = "fhPtDeltaMg";
323  dmg_name += (s+1);
324  fhPtDeltaMg[s] = new TH2F(dmg_name,dmg_name,15,0.,150.,400,-0.3,0.3);
325 
326  fOutput->Add( fhPtDeltaZg[s] );
327  fOutput->Add( fhPtDeltaRg[s] );
328  fOutput->Add( fhPtDeltaMg[s] );
329 
330  dzg_name += "BackSub";
331  fhPtDeltaZgBackSub[s] = new TH2F(dzg_name,dzg_name,15,0.,150.,200,-0.2,0.2);
332  drg_name += "BackSub";
333  fhPtDeltaRgBackSub[s] = new TH2F(drg_name,drg_name,15,0.,150.,200,r_l,r_r);
334  dmg_name += "BackSub";
335  fhPtDeltaMgBackSub[s] = new TH2F(dmg_name,dmg_name,15,0.,150.,400,-0.3,0.3);
336 
337  fOutput->Add( fhPtDeltaZgBackSub[s] );
338  fOutput->Add( fhPtDeltaRgBackSub[s] );
339  fOutput->Add( fhPtDeltaMgBackSub[s] );
340 
341  dzg_name = "fhPtDeltaZgDet";
342  fhPtDeltaZgDet[s] = new TH2F(dzg_name,dzg_name,15,0.,150.,200,-0.2,0.2);
343  drg_name = "fhPtDeltaRgDet";
344  fhPtDeltaRgDet[s] = new TH2F(drg_name,drg_name,15,0.,150.,200,r_l,r_r);
345  dmg_name = "fhPtDeltaMgDet";
346  fhPtDeltaMgDet[s] = new TH2F(dmg_name,dmg_name,15,0.,150.,400,-0.3,0.3);
347 
348  fOutput->Add( fhPtDeltaZgDet[s] );
349  fOutput->Add( fhPtDeltaRgDet[s] );
350  fOutput->Add( fhPtDeltaMgDet[s] );
351 
352  }
353 
354  fhRho = new TH1F("fhRho","fhRho",1000,0,1000);
355  fOutput->Add(fhRho);
356 
357  fhRhoSparse = new TH1F("fhRhoSparse","fhRhoSparse",1000,0,1000);
358  fOutput->Add(fhRhoSparse);
359 
360  fhPtDeltaPtAreaBackSub = new TH2F("fhPtDeltaPtAreaBackSub","fhPtDeltaPtAreaBackSub",20,0.,200.,200,-20.,20.);
362 
363  fhPtDeltaPtAreaBackSubSparse = new TH2F("fhPtDeltaPtAreaBackSubSparse","fhPtDeltaPtAreaBackSubSparse",20,0.,200.,200,-20.,20.);
365 
366  fTreeDL = new TNtuple("JetTrackTreeDL", "jet-track tree dl", "pt:eta:phi:jetm");
367  PostData(2, fTreeDL);
368 
369  fTreeDLUEBS = new TNtuple("JetTrackTreeDLUEBS", "jet-track tree dl+ue backgr sub", "pt:eta:phi:jetm");
370  PostData(3, fTreeDLUEBS);
371 
372  fTreePL = new TNtuple("JetTrackTreePL", "jet-track tree pl", "pt:eta:phi:jetm");
373  PostData(4, fTreePL);
374 
375  fRandom = new TRandom;
376 
377  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
378 }
379 
380 //________________________________________________________________________
382 
383  auto const area_nominal = TMath::Pi() * 0.4 * 0.4;
384  bool isDumpEventToTree = fRandom->Uniform() < fFractionEventsDumpedToTree;
385 
386  fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, 0.4, fastjet::E_scheme);
387  fastjet::AreaDefinition area_def(fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(0.9, 1));
388  fastjet::Selector sel_jets = fastjet::SelectorAbsEtaMax(0.9 - 0.4); //max_eta_jet
389 
390  // Fill histograms.
391 
392  std::vector<mjet> mjet_cont_pl; //probe PL
393  std::vector<mjet> mjet_cont_dlue; //response DL+UE
394  FillMjetContainer(fJetsCont1, mjet_cont_pl);
395  FillMjetContainer(fJetsCont2, mjet_cont_dlue);
396 
397 // if (fJetsCont1) {
398 // for (auto jet : fJetsCont1->accepted()) {
399 // }
400 // }
401 
402  //dump pl jets
403  if (isDumpEventToTree) {
405  PostData(4, fTreePL);
406  }
407 
408  std::vector<fastjet::PseudoJet> event_dl;
409  AddTracksToEvent(fTracksCont2, event_dl); //only pythia det level tracks
410  fastjet::ClusterSequence cs_dl(event_dl, jet_def);
411  std::vector<fastjet::PseudoJet> jets_dl = sorted_by_pt(sel_jets(cs_dl.inclusive_jets()));
412 
413  std::vector<fastjet::PseudoJet> jets_dl_filtered;
414  for (auto jet : jets_dl) {
415  if ( jet.pt() < 1. ) continue;
416  if ( jet.has_area() ) {
417  auto jarea = jet.area_4vector().perp();
418  if (jarea < (0.6 * area_nominal)) continue;
419  }
420  jets_dl_filtered.push_back(jet);
421  }
422 
423  std::vector<mjet> mjet_cont_dl;
424  FillMjetContainer(jets_dl_filtered, mjet_cont_dl);
425 
426  //PL-DL, no UE
427  std::vector<AliEmcalJet*> jets_pl_matched;
428  std::vector<fastjet::PseudoJet> jets_dl_matched;
429  for (int i = 0; i < mjet_cont_pl.size(); i++) {
430  for (int j = 0; j < mjet_cont_dl.size(); j++) {
431  auto mjet1 = mjet_cont_pl[i];
432  auto mjet2 = mjet_cont_dl[j];
433  auto dist = CalcDist(mjet1, mjet2);
434  if (dist > 0.8) continue; //loose cut
435  auto eshare = CalcEnergyShare(mjet1, mjet2);
436  if (eshare < 0.1) continue; //loose cut
437 
438  //delayed calculation of the splits
439  //deep declustering is performed
440  //only for the (roughly) matched jets
441  if (mjet1.pointerAJet) mjet1.splits = ReclusterFindHardSplits( mjet1.pointerAJet );
442  if (mjet2.pointerPJet) mjet2.splits = ReclusterFindHardSplits( *(mjet2.pointerPJet) );
443 
444  auto iscl = IsClosestPair(i,j,mjet_cont_pl,mjet_cont_dl);
445  FillResponseFromMjets(mjet1, mjet2, fhResponseDet, dist, eshare, iscl);
446  if ( (dist < 0.4) && (eshare > 0.5) && iscl ) { //strict cuts
448  if (mjet1.pointerAJet) jets_pl_matched.push_back(mjet1.pointerAJet);
449  if (mjet2.pointerPJet) jets_dl_matched.push_back( *(mjet2.pointerPJet) );
450  }
451  }
452  }
453 
454  //dump only matched jets
455  if (isDumpEventToTree) {
456  FillRespTree(jets_pl_matched, jets_dl_matched, fTreeDL);
457  PostData(2, fTreeDL);
458  }
459 
460  //FULL EVENT
461  std::vector<fastjet::PseudoJet> event_full;
462  AddTracksToEvent(fTracksCont1, event_full);
463  AddTracksToEvent(fTracksCont2, event_full);
464 
465  //get backgr-subtracted jets
466  Double_t rho;
467  Double_t rho_sparse;
468  InitializeSubtractor(event_full, rho, rho_sparse, fbcoption);
469  fhRho->Fill(rho);
470  fhRhoSparse->Fill(rho_sparse);
471 
472  std::vector<fastjet::PseudoJet> jets_backsub = GetBackSubJets(event_full);
473 
474  std::vector<fastjet::PseudoJet> jets_backsub_filtered;
475  for (auto jet : jets_backsub) {
476  if ( jet.pt() < 1. ) continue;
477  if ( jet.has_area() ) {
478  auto jarea = jet.area_4vector().perp();
479  if (jarea < (0.6 * area_nominal)) continue;
480  }
481  jets_backsub_filtered.push_back(jet);
482  }
483 
484  std::vector<mjet> mjet_container_backsub;
485  FillMjetContainer(jets_backsub_filtered, mjet_container_backsub);
486 
487  //background-subtracted jets
488  jets_pl_matched.clear();
489  std::vector<fastjet::PseudoJet> jets_backsub_filtered_matched;
490  for (int i = 0; i < mjet_cont_pl.size(); i++) {
491  for (int j = 0; j < mjet_container_backsub.size(); j++) {
492  auto mjet1 = mjet_cont_pl[i];
493  auto mjet2 = mjet_container_backsub[j];
494  auto dist = CalcDist(mjet1, mjet2);
495  if (dist > 0.8) continue; //loose cut
496  auto eshare = CalcEnergyShare(mjet1, mjet2);
497  if (eshare < 0.1) continue; //loose cut
498 
499  //delayed calculation
500  if (mjet1.pointerAJet) mjet1.splits = ReclusterFindHardSplits( mjet1.pointerAJet );
501  if (mjet2.pointerPJet) mjet2.splits = ReclusterFindHardSplits( *(mjet2.pointerPJet) );
502 
503  auto iscl = IsClosestPair(i,j,mjet_cont_pl,mjet_container_backsub);
504  FillResponseFromMjets(mjet1, mjet2, fhResponseBackSub, dist, eshare, iscl);
505  if ( (dist < 0.4) && (eshare > 0.5) && iscl ) { //strict cuts
507  if (mjet1.pointerAJet) jets_pl_matched.push_back(mjet1.pointerAJet);
508  if (mjet2.pointerPJet) jets_backsub_filtered_matched.push_back( *(mjet2.pointerPJet) );
509  }
510  }
511  }
512 
513  //and now dump only matched jets
514  if (isDumpEventToTree) {
515  FillRespTree(jets_pl_matched, jets_backsub_filtered_matched, fTreeDLUEBS);
516  PostData(3, fTreeDLUEBS);
517  }
518 
519  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
520 
521 
522  if (fCSubtractor) delete fCSubtractor;
523  if (fCSubtractorCS) delete fCSubtractorCS;
524 
525  fCSubtractor = 0;
526  fCSubtractorCS = 0;
527 
528  return kTRUE;
529 }
530 
531 //________________________________________________________________________
533 
534 
536 
539 // fTracksCont = GetParticleContainer(0); //old
540 
541 // fTracksCont1 = GetParticleContainer(0);
542 // fTracksCont1->SetClassName("AliAODMCParticle");
543 
544 // fTracksCont0 = GetParticleContainer(0);
545 // fTracksCont0->SetClassName("AliAODMCParticle");
546 
548  fTracksCont1->SetClassName("AliVTrack");
549 
551  if (fTracksCont2) fTracksCont2->SetClassName("AliVTrack");
552 
553  if (fJetsCont1 && fJetsCont1->GetArray() == 0) fJetsCont1 = 0;
554  if (fJetsCont2 && fJetsCont2->GetArray() == 0) fJetsCont2 = 0;
555 // if (fTracksCont0 && fTracksCont0->GetArray() == 0) fTracksCont0 = 0;
556  if (fTracksCont1 && fTracksCont1->GetArray() == 0) fTracksCont1 = 0;
557  if (fTracksCont2 && fTracksCont2->GetArray() == 0) fTracksCont2 = 0;
558 
559 }
560 
561 int AliAnalysisTaskSDKLResponse::FillResponseFromMjets(mjet const & mjet1, mjet const & mjet2, THnSparse* hr[4], Float_t const dist, Float_t const eshare, Bool_t const iscl) {
562 
563  std::vector<split> const & splits_1 = mjet1.splits;
564  std::vector<split> const & splits_2 = mjet2.splits;
565  auto pt1 = mjet1.pt;
566  auto pt2 = mjet2.pt;
567 
568  auto nsd_1 = splits_1.size();
569  float zg_1[4] = {0.};
570  float rg_1[4] = {0.};
571  float mg_1[4] = {-1.};
572  for (int i = 0; (i < 4) && (i < nsd_1); i++) {
573  zg_1[i] = splits_1[i].z;
574  rg_1[i] = splits_1[i].r;
575  mg_1[i] = splits_1[i].m;
576  }
577 
578  auto nsd_2 = splits_2.size();
579  float zg_2[4] = {0.};
580  float rg_2[4] = {0.};
581  float mg_2[4] = {-1.};
582  for (int i = 0; (i < 4) && (i < nsd_2); i++) {
583  zg_2[i] = splits_2[i].z;
584  rg_2[i] = splits_2[i].r;
585  mg_2[i] = splits_2[i].m;
586  }
587 
588  //dist eshare iscl pt1 zg1 rg1 mg1/pt1 nsd1 pt2 zg2 rg2 mg2/pt2
589  // 0 1 2 3 4 5 6 7 8 9 10
590 // Int_t bins[nbins] = { 10, 10, 2, 200, 20, 20, 20, 200, 20, 20, 20};
591 // Double_t xmin[nbins] = {0.0, 0.0, 0, 0.0, 0.0, 0.0, 0.0, 0., 0.0, 0.0, 0.0};
592 // Double_t xmax[nbins] = {0.8, 1.0, 2, 200., 0.5, 0.4, 0.3, 200., 0.5, 0.4, 0.3};
593 
594  Double_t xvalue[4][11]; //4 hard splits
595  for (int s = 0; s < 4; s++) {
596  xvalue[s][0] = dist;
597  xvalue[s][1] = eshare;
598  xvalue[s][2] = 0.1 + iscl;
599 
600  xvalue[s][3] = pt1;
601  xvalue[s][4] = zg_1[s];
602  xvalue[s][5] = rg_1[s];
603  xvalue[s][6] = mg_1[s] / pt1;
604 
605  xvalue[s][7] = pt2;
606  xvalue[s][8] = zg_2[s];
607  xvalue[s][9] = rg_2[s];
608  xvalue[s][10] = mg_2[s] / pt2;
609 
610  hr[s]->Fill( xvalue[s] );
611  }
612 
613  return 0;
614 
615 }
616 
617 
618 int AliAnalysisTaskSDKLResponse::FillDeltasFromMjets(mjet const & mjet1, mjet const & mjet2, TH2F *hptdpt, TH2F *h1[4], TH2F *h2[4], TH2F *h3[4]) {
619 
620  auto dist = CalcDist(mjet1, mjet2);
621  if ( dist > (0.4*0.8) ) return 1; //tight cut - only "matched" jets
622 
623  auto eshare = CalcEnergyShare(mjet1, mjet2);
624  if (eshare < 0.5) return 1; //eshare cut
625 
626  std::vector<split> const & splits_1 = mjet1.splits;
627  std::vector<split> const & splits_2 = mjet2.splits;
628  auto pt1 = mjet1.pt;
629  auto pt2 = mjet2.pt;
630 
631  int nsd_1 = splits_1.size();
632  float zg_1[4] = {0.};
633  float rg_1[4] = {0.};
634  float mg_1[4] = {-1.};
635  for (int i = 0; (i < 4) && (i < nsd_1); i++) {
636  zg_1[i] = splits_1[i].z;
637  rg_1[i] = splits_1[i].r;
638  mg_1[i] = splits_1[i].m;
639  }
640 
641  int nsd_2 = splits_2.size();
642  float zg_2[4] = {0.};
643  float rg_2[4] = {0.};
644  float mg_2[4] = {-1.};
645  for (int i = 0; (i < 4) && (i < nsd_2); i++) {
646  zg_2[i] = splits_2[i].z;
647  rg_2[i] = splits_2[i].r;
648  mg_2[i] = splits_2[i].m;
649  }
650 
651  hptdpt->Fill(pt1, pt2 - pt1);
652  for (int s = 0; s < 4; s++) {
653  if ( (zg_1[s] > 0.05) && (zg_2[s] > 0.05) ) {
654  h1[s]->Fill(pt1, zg_2[s] - zg_1[s]);
655  h2[s]->Fill(pt1, rg_2[s] - rg_1[s]);
656  h3[s]->Fill(pt1, mg_2[s] - mg_1[s]);
657  }//only tagged jets
658  }//splits iter
659 
660  return 0;
661 
662 }
663 
664 void AliAnalysisTaskSDKLResponse::FillMjetContainer(AliJetContainer *jet_container, std::vector <mjet> & mjet_container) {
665  if (jet_container) {
666  for ( auto jet : jet_container->accepted() ) {
667 
668  double pt_scalar = 0.0;
669  std::vector<double> track_pts;
670  std::vector<int> track_labels;
671 
672  UShort_t ntracks = jet->GetNumberOfTracks();
673  for (int j = 0; j < ntracks; j++) {
674  auto jtrack = jet->Track(j);
675  auto track_pt = jtrack->Pt();
676  pt_scalar += track_pt;
677  auto mclabel = jtrack->GetLabel();
678  if (mclabel > -1) { //keep only MC labels
679  track_pts.push_back( track_pt );
680  track_labels.push_back( mclabel );
681  }
682  }
683  std::vector<split> empty;
684  mjet_container.push_back( mjet{ jet->Pt(), pt_scalar, jet->Eta(), jet->Phi(), jet->Area(), track_pts, track_labels, empty, jet, nullptr } );
685  }
686  }
687 }
688 
689 void AliAnalysisTaskSDKLResponse::FillMjetContainer(std::vector <fastjet::PseudoJet> & jet_container, std::vector <mjet> & mjet_container) {
690 
691  for ( int k = 0; k < jet_container.size(); k++ ) {
692  auto jet = jet_container[k];
693  if ( jet.perp() < 1.e-6 ) continue; //skip ghost-only jets
694  double area = 0.0;
695  double pt_scalar = 0.0;
696  std::vector<double> track_pts;
697  std::vector<int> track_labels;
698 
699  for (auto c : jet.constituents() ) {
700  pt_scalar += c.pt();
701  auto mclabel = c.user_index();
702  if (mclabel > -1) { //keep only MC labels
703  track_pts.push_back( c.pt() );
704  track_labels.push_back( mclabel );
705  }
706  }
707 
708  if ( jet.has_area() ) area = jet.area_4vector().perp();
709  std::vector<split> empty;
710  mjet_container.push_back( mjet{ jet.perp(), pt_scalar, jet.eta(), jet.phi(), area, track_pts, track_labels, empty, nullptr, &(jet_container[k]) } );
711  }
712 
713 }
714 
715 //________________________________________________________________________
717 {
718  // Run analysis code here, if needed. It will be executed before FillHistograms().
719 
720  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
721 }
722 
723 //________________________________________________________________________
725 {
726  // Called once at the end of the analysis.
727 }
728 
730 
731  auto phi1 = mjet1.phi;
732  auto eta1 = mjet1.eta;
733  auto phi2 = mjet2.phi;
734  auto eta2 = mjet2.eta;
735 
736  auto phi_diff = fabs(phi1 - phi2);
737  if ( phi_diff > TMath::Pi() ) phi_diff = 2.0 * TMath::Pi() - phi_diff;
738  auto eta_diff = fabs(eta1 - eta2);
739 
740  auto dist = sqrt(phi_diff*phi_diff + eta_diff*eta_diff);
741 
742  return dist;
743 
744 }
745 
746 Bool_t AliAnalysisTaskSDKLResponse::IsClosestPair(int const idx1, int const idx2, std::vector <mjet> const & mjet_container1, std::vector <mjet> const & mjet_container2) {
747 
748  double min_dist12 = 1000.0;
749  int min_dist12_index2 = -7;
750 
751  auto size2 = mjet_container2.size();
752  for (int k = 0; k < size2; k++) {
753  auto cur_dist = CalcDist( mjet_container1[idx1], mjet_container2[k] );
754  if (cur_dist < min_dist12) {
755  min_dist12 = cur_dist;
756  min_dist12_index2 = k;
757  }
758  }
759 
760  if ( idx2 != min_dist12_index2 ) return kFALSE;
761 
762  double min_dist21 = 1000.0;
763  int min_dist21_index1 = -7;
764 
765  auto size1 = mjet_container1.size();
766  for (int l = 0; l < size1; l++) {
767  auto cur_dist = CalcDist( mjet_container2[idx2], mjet_container1[l] );
768  if (cur_dist < min_dist21) {
769  min_dist21 = cur_dist;
770  min_dist21_index1 = l;
771  }
772  }
773 
774  if ( idx1 != min_dist21_index1 ) return kFALSE;
775 
776  return kTRUE;
777 
778 }
779 
781 
782  auto track_pts2 = mjet2.track_pts;
783  auto track_labels1 = mjet1.track_labels;
784  auto track_labels2 = mjet2.track_labels;
785 
786  double pt_tot_from_jet1 = 0.0;
787  for (int j = 0; j < track_labels2.size(); j++) {
788  if ( std::find( track_labels1.begin(), track_labels1.end(), track_labels2[j] ) != track_labels1.end() ) {
789  pt_tot_from_jet1 += track_pts2[j];
790  }
791  }
792 
793  return pt_tot_from_jet1/mjet2.pt_scalar;
794 
795 }
796 
797 void AliAnalysisTaskSDKLResponse::FillRespTree(std::vector<AliEmcalJet*> const & probe_jets, std::vector<fastjet::PseudoJet> const & resp_jets, TNtuple* tree) {
798 
799  int njets = probe_jets.size();
800 
801  for (int i = 0; i < njets; i++) {
802 
803  //probe
804  auto pjet = probe_jets[i];
805  if ( pjet->Pt() < 10.) continue; //cut only on the probe jet
806 
807  UShort_t ntracks = pjet->GetNumberOfTracks();
808  tree->Fill(pjet->Pt(), pjet->Eta(), pjet->Phi(), ntracks);
809  for (int j = 0; j < ntracks; j++) {
810  auto jtrack = pjet->Track(j);
811  if (jtrack->Pt() > 1.e-5) {
812  tree->Fill(jtrack->Pt(), jtrack->Eta(), jtrack->Phi(), -108);
813  }
814  }
815 
816  //response
817  auto rjet = resp_jets[i];
818  int nconst = 0;
819  for (auto c : rjet.constituents() ) {
820  if ( c.pt() > 1.e-5 ) nconst++;
821  }
822  tree->Fill(rjet.pt(), rjet.eta(), rjet.phi(), nconst);
823  for ( auto c : rjet.constituents() ) {
824  if ( c.pt() > 1.e-5 ) {
825  tree->Fill(c.pt(), c.eta(), c.phi(), -7);
826  }
827  }
828 
829  }
830 
831 }
THnSparse * fhResponseBackSub[4]
! distribution of all
Bool_t IsClosestPair(int const i, int const j, std::vector< mjet > const &mjet_container1, std::vector< mjet > const &mjet_container2)
int FillResponseFromMjets(mjet const &mjet1, mjet const &mjet2, THnSparse *hr[4], Float_t const dist, Float_t const eshare, Bool_t const iscl)
double Double_t
Definition: External.C:58
Definition: External.C:236
AliJetContainer * GetJetContainer(Int_t i=0) const
void SetLeadingHadronType(Int_t t)
TCanvas * c
Definition: TestFitELoss.C:172
std::vector< double > track_pts
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
THnSparse * fhResponse[4]
! distribution of all
void FillRespTree(std::vector< AliEmcalJet * > const &probe_jets, std::vector< fastjet::PseudoJet > const &resp_jets, TNtuple *tree)
void SetPercAreaCut(Float_t p)
int FillDeltasFromMjets(mjet const &mjet1, mjet const &mjet2, TH2F *hptdpt, TH2F *h1[4], TH2F *h2[4], TH2F *h3[4])
Int_t nCentBins
Float_t CalcEnergyShare(mjet const &mjet1, mjet const &mjet2)
fastjet::ClusterSequenceArea * fCSubtractorCS
Float_t CalcDist(mjet const &mjet1, mjet const &mjet2)
Container for particles within the EMCAL framework.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
void SetRhoName(const char *n)
THnSparse * fhResponseDet[4]
! distribution of all
void AddTracksToEvent(AliParticleContainer *cont, std::vector< fastjet::PseudoJet > &event)
int Int_t
Definition: External.C:63
void SetJetPtCut(Float_t cut)
float Float_t
Definition: External.C:68
int InitializeSubtractor(std::vector< fastjet::PseudoJet > const &event_full, Double_t &rho, Double_t &rho_sparse, Int_t opt=0)
Bool_t FillHistograms()
Function filling histograms.
TClonesArray * GetArray() const
std::vector< split > ReclusterFindHardSplits(AliEmcalJet *jet)
void ExecOnce()
Perform steps needed to initialize the analysis.
AliParticleContainer * fTracksCont2
Tracks.
void SetClassName(const char *clname)
std::vector< split > splits
AliParticleContainer * fTracksCont1
Jets.
void FillTree(std::vector< fastjet::PseudoJet > const &jets, TNtuple *tree)
AliEmcalList * fOutput
!output list
std::vector< int > track_labels
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
static AliAnalysisTaskSDKLResponse * AddTaskSoftDropResponse(const char *ntracks="usedefault", const char *njets1="Jets1", const char *njets2="Jets2", const char *nrho="Rho", Int_t nCentBins=1, Double_t jetradius=0.4, Double_t jetptcut=1, Double_t jetareacut=0.6, const char *type="EMCAL", Int_t backgroption=0, Int_t leadhadtype=0, Double_t fractioneventsfortree=1.e-6, const char *taskname="AliAnalysisTaskSDKLResponse")
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
std::vector< fastjet::PseudoJet > GetBackSubJets(std::vector< fastjet::PseudoJet > const &event_full)
TTree * tree
void FillMjetContainer(AliJetContainer *jet_container, std::vector< mjet > &mjet_container)
Container for jet within the EMCAL jet framework.
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
fastjet::contrib::ConstituentSubtractor * fCSubtractor