AliPhysics  32e057f (32e057f)
AliAnalysisTaskJetCoreEmcal.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, 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 #include <TClonesArray.h>
17 #include <TList.h>
18 #include "TChain.h"
19 #include "TTree.h"
20 #include "TMath.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 #include "TH3F.h"
24 #include "THnSparse.h"
25 #include "TCanvas.h"
26 #include "TRandom3.h"
27 #include "AliLog.h"
28 
29 #include <AliAnalysisManager.h>
30 #include <AliVEventHandler.h>
31 #include <AliVEvent.h>
32 #include <AliVCluster.h>
33 #include <AliVParticle.h>
34 #include <AliLog.h>
35 
36 #include "AliTLorentzVector.h"
37 #include "AliEmcalJet.h"
38 #include "AliRhoParameter.h"
39 #include "AliJetContainer.h"
40 #include "AliParticleContainer.h"
41 #include "AliClusterContainer.h"
42 
44 
48 
54  fHistManager(),
55  fCentMin(0.),
56  fCentMax(100.),
57  fTTLowRef(8.),
58  fTTUpRef(9.),
59  fTTLowSig(20.),
60  fTTUpSig(50.),
61  fNRPBins(50),
62  fFrac(0.8),
63  fJetEtaMin(-.5),
64  fJetEtaMax(.5),
65  fJetHadronDeltaPhi(0.6),
66  fJetContName(""),
67  fRandom(0),
68  fHistEvtSelection(0x0),
69  fHJetSpec(0x0),
70  fh1TrigRef(0x0),
71  fh1TrigSig(0x0),
72  fh2Ntriggers(0x0),
73  fh2RPJetsC10(0x0),
74  fh2RPJetsC20(0x0),
75  fh2RPTC10(0x0),
76  fh2RPTC20(0x0)
77 
78 {
79 }
80 
87  AliAnalysisTaskEmcalJet(name, kTRUE),
88  fHistManager(name),
89  fCentMin(0.),
90  fCentMax(100.),
91  fTTLowRef(8),
92  fTTUpRef(9.),
93  fTTLowSig(20.),
94  fTTUpSig(50.),
95  fNRPBins(50),
96  fFrac(0.8),
97  fJetEtaMin(-.5),
98  fJetEtaMax(.5),
99  fJetHadronDeltaPhi(0.6),
100  fJetContName(""),
101  fRandom(0),
102  fHistEvtSelection(0x0),
103  fHJetSpec(0x0),
104  fh1TrigRef(0x0),
105  fh1TrigSig(0x0),
106  fh2Ntriggers(0x0),
107  fh2RPJetsC10(0x0),
108  fh2RPJetsC20(0x0),
109  fh2RPTC10(0x0),
110  fh2RPTC20(0x0)
111 {
113 }
114 
119 {
120 }
121 
127 {
129 
135 
136  TIter next(fHistManager.GetListOfHistograms());
137  TObject* obj = 0;
138  while ((obj = next())) {
139  fOutput->Add(obj);
140  }
141 
142  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
143 }
144 
145 /*
146  * This function allocates the histograms for basic EMCal cluster QA.
147  * A set of histograms (energy, eta, phi, number of cluster) is allocated
148  * per each cluster container and per each centrality bin.
149  */
151 {
152  TString histname;
153  TString histtitle;
154  TString groupname;
155  AliClusterContainer* clusCont = 0;
156  TIter next(&fClusterCollArray);
157  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
158  groupname = clusCont->GetName();
159  // Protect against creating the histograms twice
160  if (fHistManager.FindObject(groupname)) {
161  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
162  continue;
163  }
164  fHistManager.CreateHistoGroup(groupname);
165  for (Int_t cent = 0; cent < fNcentBins; cent++) {
166  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
167  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
168  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
169 
170  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
171  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
172  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
173 
174  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
175  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
176  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
177 
178  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
179  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
180  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
181 
182  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
183  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
184  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
185 
186  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
187  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
188  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
189 
190  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
191  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
192  if (fForceBeamType != kpp) {
193  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
194  }
195  else {
196  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
197  }
198  }
199  }
200 
201  histname = "fHistSumNClusters";
202  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
203  if (fForceBeamType != kpp) {
204  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
205  }
206  else {
207  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
208  }
209 }
210 
211 /*
212  * This function allocates the histograms for basic EMCal QA.
213  * One 2D histogram with the cell energy spectra and the number of cells
214  * per event is allocated per each centrality bin.
215  */
217 {
218  TString histname;
219  TString histtitle;
220  TString groupname(fCaloCellsName);
221 
222  fHistManager.CreateHistoGroup(groupname);
223  for (Int_t cent = 0; cent < fNcentBins; cent++) {
224  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
225  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
226  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
227 
228  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
229  histtitle = TString::Format("%s;number of cells;events", histname.Data());
230  if (fForceBeamType != kpp) {
231  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
232  }
233  else {
234  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
235  }
236  }
237 }
238 
239 /*
240  * This function allocates the histograms for basic tracking QA.
241  * A set of histograms (pT, eta, phi, difference between kinematic properties
242  * at the vertex and at the EMCal surface, number of tracks) is allocated
243  * per each particle container and per each centrality bin.
244  */
246 {
247  TString histname;
248  TString histtitle;
249  TString groupname;
250  AliParticleContainer* partCont = 0;
251  TIter next(&fParticleCollArray);
252  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
253  groupname = partCont->GetName();
254  // Protect against creating the histograms twice
255  if (fHistManager.FindObject(groupname)) {
256  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
257  continue;
258  }
259  fHistManager.CreateHistoGroup(groupname);
260  for (Int_t cent = 0; cent < fNcentBins; cent++) {
261  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
262  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
263  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
264 
265  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
266  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
267  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
268 
269  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
270  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
271  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
272 
273  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
274  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
275  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
276  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
277 
278  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
279  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
280  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
281 
282  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
283  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{p}_{T,track}^{vertex} - #it{p}_{T,track}^{EMCal} (GeV/#it{c});counts", histname.Data());
284  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
285 
286  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
287  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
288  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
289  }
290 
291  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
292  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
293  if (fForceBeamType != kpp) {
294  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
295  }
296  else {
297  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
298  }
299  }
300  }
301 
302  histname = "fHistSumNTracks";
303  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
304  if (fForceBeamType != kpp) {
305  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
306  }
307  else {
308  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
309  }
310 }
311 
312 /*
313  * This function allocates the histograms for basic jet QA.
314  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
315  * per each jet container and per each centrality bin.
316  */
318 {
319  TString histname;
320  TString histtitle;
321  TString groupname;
322  AliJetContainer* jetCont = 0;
323  TIter next(&fJetCollArray);
324  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
325  groupname = jetCont->GetName();
326  // Protect against creating the histograms twice
327  if (fHistManager.FindObject(groupname)) {
328  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
329  continue;
330  }
331  fHistManager.CreateHistoGroup(groupname);
332  for (Int_t cent = 0; cent < fNcentBins; cent++) {
333  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
334  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
335  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
336 
337  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
338  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
339  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
340 
341  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
342  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
343  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
344 
345  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
346  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
347  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
348 
349  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
350  histtitle = TString::Format("%s;number of jets;events", histname.Data());
351  if (fForceBeamType != kpp) {
352  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
353  }
354  else {
355  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
356  }
357 
358  if (!jetCont->GetRhoName().IsNull()) {
359  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
360  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
361  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
362  }
363  }
364  }
365 }
366 
368 {
369 
370  Bool_t oldStatus = TH1::AddDirectoryStatus();
371  TH1::AddDirectory(kFALSE);
372 
373  // set seed
374  fRandom = new TRandom3(0);
375 
376  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
377  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
378  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
379  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
380  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
381  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
382  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
383 
384  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
385  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
386  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
387  fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
388  fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
389  fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
390  fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
391 
393 
394  fOutput->Add(fh1TrigRef);
395  fOutput->Add(fh1TrigSig);
396  fOutput->Add(fh2Ntriggers);
397  fOutput->Add(fh2RPJetsC10);
398  fOutput->Add(fh2RPJetsC20);
399  fOutput->Add(fh2RPTC10);
400  fOutput->Add(fh2RPTC20);
401 
402  const Int_t dimSpec = 5;
403  const Int_t nBinsSpec[dimSpec] = {100,100, 140, 50, fNRPBins};
404  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
405  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50, static_cast<Double_t>(fNRPBins)};
406  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
407 
408  // comment out since I want finer binning in jet area, to make it easier
409  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
410  // for inclusive jets)
411  //change binning in jet area
412 // Double_t *xPt6 = new Double_t[7];
413 // xPt6[0] = 0.;
414 // xPt6[1]=0.07;
415 // xPt6[2]=0.2;
416 // xPt6[3]=0.4;
417 // xPt6[4]=0.6;
418 // xPt6[5]=0.8;
419 // xPt6[6]=1;
420 // fHJetSpec->SetBinEdges(1,xPt6);
421 // delete [] xPt6;
422 
423  fOutput->Add(fHJetSpec);
424 
425  // =========== Switch on Sumw2 for all histos ===========
426  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
427  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
428  if (h1){
429  h1->Sumw2();
430  continue;
431  }
432  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
433  if (hn){
434  hn->Sumw2();
435  }
436  }
437 
438  // add QA plots from fEventCuts
439  fEventCuts.AddQAplotsToList(fOutput);
440 
441  TH1::AddDirectory(oldStatus);
442 
443  PostData(1, fOutput);
444 }
445 
453 {
454 
455  fHistEvtSelection->Fill(1); // number of events before event selection
456  AliVEvent *ev = InputEvent();
457  if (!fEventCuts.AcceptEvent(ev)) {
458  fHistEvtSelection->Fill(2);
459  return kTRUE;
460  }
461 
462  DoJetLoop();
463  DoTrackLoop();
464  //DoClusterLoop();
465  //DoCellLoop();
466  DoJetCoreLoop();
467 
468  return kTRUE;
469 }
470 
472 {
473 
474  // Do jet core analysis and fill histograms.
476  if(!jetCont) {
477  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
478  TIter next(&fJetCollArray);
479  while ((jetCont = static_cast<AliJetContainer*>(next())))
480  AliError(Form("%s",jetCont->GetName()));
481  AliFatal("Exit...");
482  return;
483  }
484 
485  // centrality selection
486  if(fDebug) Printf("centrality: %f\n", fCent);
487  if ((fCent>fCentMax) || (fCent<fCentMin)) {
488  fHistEvtSelection->Fill(4);
489  return;
490  }
491  fHistEvtSelection->Fill(0);
492 
493  // Background
494  Double_t rho = 0;
495  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
496  if(fDebug) Printf("rho = %f, rho check = %f",rho, GetRhoVal(0));
497 
498  // Choose trigger track
499  Int_t nT=0;
500  TList ParticleList;
501  Double_t minT=0;
502  Double_t maxT=0;
503  Int_t number=0;
504  Double_t dice=fRandom->Uniform(0,1);
505  if(dice>fFrac){
506  minT=fTTLowRef;
507  maxT=fTTUpRef;
508  }
509  if(dice<=fFrac){
510  minT=fTTLowSig;
511  maxT=fTTUpSig;
512  }
513  nT=SelectTrigger(&ParticleList,minT,maxT,number);
514  if(fDebug) Printf("%s class ---> n triggers between %f and %f = %i, index of trigger chosen = %i",dice>fFrac?"ref.":"sig.",minT,maxT,number,nT);
515  if(nT<0) return;
516 
517  if(dice>fFrac) fh1TrigRef->Fill(number);
518  if(dice<=fFrac)fh1TrigSig->Fill(number);
519 
520 
521  // particle loop -
522  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
523  // histogram
524  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
525  if(tt!=nT) continue;
526  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
527  if(!partback) continue;
528  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
529  // if(partback->Pt()<8) continue;
530 
531  Int_t injet4=0;
532  Int_t injet=0;
533 
534  fh2Ntriggers->Fill(fCent,partback->Pt());
535  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
536  if(fCent<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
537  if(fCent<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
538 
539  Double_t etabig=0;
540  Double_t ptbig=0;
541  Double_t areabig=0;
542  Double_t phibig=0.;
543  // Double_t areasmall=0;
544 
545  TString histname;
546  TString groupname;
547  groupname = jetCont->GetName();
548  UInt_t count = 0;
549  for(auto jetbig : jetCont->accepted()) {
550  if (!jetbig) continue;
551  count++;
552  ptbig = jetbig->Pt();
553  etabig = jetbig->Eta();
554  phibig = jetbig->Phi();
555  if(ptbig==0) continue;
556  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
557  areabig = jetbig->Area();
558  Double_t ptcorr=ptbig-rho*areabig;
559  //JJJ - perhaps should change eta selection if implemented in jet container
560  if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
561  if(areabig>=0.07) injet=injet+1;
562  if(areabig>=0.4) injet4=injet4+1;
563  Double_t dphi=RelativePhi(partback->Phi(),phibig);
564  if(fDebug) Printf("jet properties...\n\teta = %f \t phi = %f \t pt = %f \t relativephi = %f\t area = %f\t rho = %f",etabig,phibig,ptbig,dphi,areabig,rho);
565  // selection on relative phi
566  if(TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
567 
568  if(fCent<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
569  if(fCent<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
570 
571  Float_t phitt=partback->Phi();
572  if(phitt<0)phitt+=TMath::Pi()*2.;
573  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
574 
575  Double_t fillspec[] = {fCent,jetbig->Area(),ptcorr,partback->Pt(), static_cast<Double_t>(phiBintt)};
576  fHJetSpec->Fill(fillspec);
577  }
578  }
579 }
580 
586 {
587  TString histname;
588  TString groupname;
589  AliJetContainer* jetCont = 0;
590  TIter next(&fJetCollArray);
591  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
592  groupname = jetCont->GetName();
593  UInt_t count = 0;
594  for(auto jet : jetCont->accepted()) {
595  if (!jet) continue;
596  count++;
597 
598  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
599  fHistManager.FillTH1(histname, jet->Pt());
600 
601  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
602  fHistManager.FillTH1(histname, jet->Area());
603 
604  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
605  fHistManager.FillTH1(histname, jet->Phi());
606 
607  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
608  fHistManager.FillTH1(histname, jet->Eta());
609 
610  if (jetCont->GetRhoParameter()) {
611  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
612  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
613  }
614  }
615  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
616  fHistManager.FillTH1(histname, count);
617  }
618 }
619 
625 {
627 
628  TString histname;
629  TString groupname;
630  UInt_t sumAcceptedTracks = 0;
631  AliParticleContainer* partCont = 0;
632  TIter next(&fParticleCollArray);
633  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
634  groupname = partCont->GetName();
635  UInt_t count = 0;
636  for(auto part : partCont->accepted()) {
637  if (!part) continue;
638  count++;
639 
640  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
641  fHistManager.FillTH1(histname, part->Pt());
642 
643  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
644  fHistManager.FillTH1(histname, part->Phi());
645 
646  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
647  fHistManager.FillTH1(histname, part->Eta());
648 
649  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
650  const AliVTrack* track = static_cast<const AliVTrack*>(part);
651 
652  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
653  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
654 
655  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
656  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
657 
658  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
659  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
660 
661  if (clusCont) {
662  Int_t iCluster = track->GetEMCALcluster();
663  if (iCluster >= 0) {
664  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
665  if (cluster) {
666  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
667  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
668  }
669  }
670  }
671  }
672  }
673  sumAcceptedTracks += count;
674 
675  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
676  fHistManager.FillTH1(histname, count);
677  }
678 
679  histname = "fHistSumNTracks";
680  fHistManager.FillTH1(histname, sumAcceptedTracks);
681 }
682 
688 {
689  TString histname;
690  TString groupname;
691  UInt_t sumAcceptedClusters = 0;
692  AliClusterContainer* clusCont = 0;
693  TIter next(&fClusterCollArray);
694  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
695  groupname = clusCont->GetName();
696 
697  for(auto cluster : clusCont->all()) {
698  if (!cluster) continue;
699 
700  if (cluster->GetIsExotic()) {
701  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
702  fHistManager.FillTH1(histname, cluster->E());
703  }
704  }
705 
706  UInt_t count = 0;
707  for(auto cluster : clusCont->accepted()) {
708  if (!cluster) continue;
709  count++;
710 
711  AliTLorentzVector nPart;
712  cluster->GetMomentum(nPart, fVertex);
713 
714  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
715  fHistManager.FillTH1(histname, cluster->E());
716 
717  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
718  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
719 
720  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
721  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
722 
723  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
724  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
725 
726  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
727  fHistManager.FillTH1(histname, nPart.Eta());
728  }
729  sumAcceptedClusters += count;
730 
731  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
732  fHistManager.FillTH1(histname, count);
733  }
734 
735  histname = "fHistSumNClusters";
736  fHistManager.FillTH1(histname, sumAcceptedClusters);
737 }
738 
744 {
745  if (!fCaloCells) return;
746 
747  TString histname;
748 
749  const Short_t ncells = fCaloCells->GetNumberOfCells();
750 
751  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
752  fHistManager.FillTH1(histname, ncells);
753 
754  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
755  for (Short_t pos = 0; pos < ncells; pos++) {
756  Double_t amp = fCaloCells->GetAmplitude(pos);
757 
758  fHistManager.FillTH1(histname, amp);
759  }
760 }
761 
767 {
769 }
770 
779 {
780  return kTRUE;
781 }
782 
787 {
788 }
789 
790 THnSparse* AliAnalysisTaskJetCoreEmcal::NewTHnSparseF(const char* name, UInt_t entries)
791 {
792  // generate new THnSparseF, axes are defined in GetDimParams()
793 
794  Int_t count = 0;
795  UInt_t tmp = entries;
796  while(tmp!=0){
797  count++;
798  tmp = tmp &~ -tmp; // clear lowest bit
799  }
800 
801  TString hnTitle(name);
802  const Int_t dim = count;
803  Int_t nbins[dim];
804  Double_t xmin[dim];
805  Double_t xmax[dim];
806 
807  Int_t i=0;
808  Int_t c=0;
809  while(c<dim && i<32){
810  if(entries&(1<<i)){
811 
812  TString label("");
813  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
814  hnTitle += Form(";%s",label.Data());
815  c++;
816  }
817 
818  i++;
819  }
820  hnTitle += ";";
821 
822  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
823 }
824 
826 {
827 
828  // stores label and binning of axis for THnSparse
829 
830  const Double_t pi = TMath::Pi();
831  switch(iEntry){
832  case 0:
833  label = "V0 centrality (%)";
834  nbins = 10;
835  xmin = 0.;
836  xmax = 100.;
837  break;
838  case 1:
839  label = "corrected jet pt";
840  nbins = 20;
841  xmin = 0.;
842  xmax = 200.;
843  break;
844  case 2:
845  label = "track pT";
846  nbins = 9;
847  xmin = 0.;
848  xmax = 150;
849  break;
850  case 3:
851  label = "deltaR";
852  nbins = 15;
853  xmin = 0.;
854  xmax = 1.5;
855  break;
856  case 4:
857  label = "deltaEta";
858  nbins = 8;
859  xmin = -1.6;
860  xmax = 1.6;
861  break;
862  case 5:
863  label = "deltaPhi";
864  nbins = 90;
865  xmin = -0.5*pi;
866  xmax = 1.5*pi;
867  break;
868  case 6:
869  label = "leading track";
870  nbins = 13;
871  xmin = 0;
872  xmax = 50;
873  break;
874  case 7:
875  label = "trigger track";
876  nbins =10;
877  xmin = 0;
878  xmax = 50;
879  break;
880  }
881 }
882 
884 
885  Int_t index=-1;
886  Int_t triggers[100];
887 
888  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
889 
890  Int_t im=0;
891 
892  TString groupname = "";
894  groupname = partCont->GetName();
895  UInt_t iCount = 0;
896  for(auto part : partCont->accepted()) {
897  if (!part) continue;
898  list->Add(part);
899  iCount++;
900  if(part->Pt()>=minT && part->Pt()<maxT){
901  triggers[im]=iCount-1;
902  im=im+1;
903  }
904  }
905  number=im;
906  Int_t rd=0;
907  if(im>0) rd=fRandom->Integer(im);
908  index=triggers[rd];
909  return index;
910 }
911 
913 
914  // get relative DeltaPhi from -pi < DeltaPhi < pi
915 
916  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
917  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
918  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
919  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
920  double dphi = mphi-vphi;
921  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
922  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
923  return dphi;//dphi in [-Pi, Pi]
924 }
925 
927 {
928  Int_t phibin=-1;
929  if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
930  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
931  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
932  if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
933  return phibin;
934 }
THashList * CreateHistoGroup(const char *groupname)
Create a new group of histograms within a parent group.
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Definition: External.C:236
Int_t SelectTrigger(TList *list, Double_t minT, Double_t maxT, Int_t &number)
AliJetContainer * GetJetContainer(Int_t i=0) const
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
TCanvas * c
Definition: TestFitELoss.C:172
Int_t fCentBin
!event centrality bin
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const AliClusterIterableContainer all() const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
float Float_t
Definition: External.C:68
Double_t Phi_0_2pi() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Definition: External.C:212
AliVCluster * GetAcceptCluster(Int_t i) const
const AliClusterIterableContainer accepted() const
THnSparse * NewTHnSparseF(const char *name, UInt_t entries)
Double_t fCent
!event centrality
TString fCaloCellsName
name of calo cell collection
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TObjArray fJetCollArray
jet collection array
Double_t RelativePhi(Double_t mphi, Double_t vphi)
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
Double_t GetRhoVal(Int_t i=0) const
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
THistManager fHistManager
Histogram manager.
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
const AliParticleIterableContainer accepted() const
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const Double_t pi
const Int_t nbins
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
void GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
Container structure for EMCAL clusters.
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Int_t fNbins
no. of pt bins