AliPhysics  2f4fa0a (2f4fa0a)
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  fFillTrackHistograms(kTRUE),
68  fFillJetHistograms(kTRUE),
69  fRandom(0),
70  fHistEvtSelection(0x0),
71  fHJetSpec(0x0),
72  fh1TrigRef(0x0),
73  fh1TrigSig(0x0),
74  fh2Ntriggers(0x0),
75  fh2RPJetsC10(0x0),
76  fh2RPJetsC20(0x0),
77  fh2RPTC10(0x0),
78  fh2RPTC20(0x0),
79  fHJetPhiCorr(0x0),
80  fhDphiPtSig(0x0),
81  fhDphiPtRef(0x0)
82 {
83 }
84 
91  AliAnalysisTaskEmcalJet(name, kTRUE),
92  fHistManager(name),
93  fCentMin(0.),
94  fCentMax(100.),
95  fTTLowRef(8),
96  fTTUpRef(9.),
97  fTTLowSig(20.),
98  fTTUpSig(50.),
99  fNRPBins(50),
100  fFrac(0.8),
101  fJetEtaMin(-.5),
102  fJetEtaMax(.5),
103  fJetHadronDeltaPhi(0.6),
104  fJetContName(""),
105  fFillTrackHistograms(kTRUE),
106  fFillJetHistograms(kTRUE),
107  fRandom(0),
108  fHistEvtSelection(0x0),
109  fHJetSpec(0x0),
110  fh1TrigRef(0x0),
111  fh1TrigSig(0x0),
112  fh2Ntriggers(0x0),
113  fh2RPJetsC10(0x0),
114  fh2RPJetsC20(0x0),
115  fh2RPTC10(0x0),
116  fh2RPTC20(0x0),
117  fHJetPhiCorr(0x0),
118  fhDphiPtSig(0x0),
119  fhDphiPtRef(0x0)
120 {
122 }
123 
128 {
129 }
130 
136 {
138 
141  //AllocateClusterHistograms();
142  //AllocateCellHistograms();
144 
145  TIter next(fHistManager.GetListOfHistograms());
146  TObject* obj = 0;
147  while ((obj = next())) {
148  fOutput->Add(obj);
149  }
150 
151  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
152 }
153 
154 /*
155  * This function allocates the histograms for basic EMCal cluster QA.
156  * A set of histograms (energy, eta, phi, number of cluster) is allocated
157  * per each cluster container and per each centrality bin.
158  */
160 {
161  TString histname;
162  TString histtitle;
163  TString groupname;
164  AliClusterContainer* clusCont = 0;
165  TIter next(&fClusterCollArray);
166  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
167  groupname = clusCont->GetName();
168  // Protect against creating the histograms twice
169  if (fHistManager.FindObject(groupname)) {
170  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
171  continue;
172  }
173  fHistManager.CreateHistoGroup(groupname);
174  for (Int_t cent = 0; cent < fNcentBins; cent++) {
175  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
176  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
177  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
178 
179  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
180  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
181  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
182 
183  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
184  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
185  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
186 
187  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
188  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
189  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
190 
191  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
192  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
193  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
194 
195  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
196  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
197  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
198 
199  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
200  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
201  if (fForceBeamType != kpp) {
202  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
203  }
204  else {
205  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
206  }
207  }
208  }
209 
210  histname = "fHistSumNClusters";
211  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
212  if (fForceBeamType != kpp) {
213  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
214  }
215  else {
216  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
217  }
218 }
219 
220 /*
221  * This function allocates the histograms for basic EMCal QA.
222  * One 2D histogram with the cell energy spectra and the number of cells
223  * per event is allocated per each centrality bin.
224  */
226 {
227  TString histname;
228  TString histtitle;
229  TString groupname(fCaloCellsName);
230 
231  fHistManager.CreateHistoGroup(groupname);
232  for (Int_t cent = 0; cent < fNcentBins; cent++) {
233  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
234  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
235  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
236 
237  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
238  histtitle = TString::Format("%s;number of cells;events", histname.Data());
239  if (fForceBeamType != kpp) {
240  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
241  }
242  else {
243  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
244  }
245  }
246 }
247 
248 /*
249  * This function allocates the histograms for basic tracking QA.
250  * A set of histograms (pT, eta, phi, difference between kinematic properties
251  * at the vertex and at the EMCal surface, number of tracks) is allocated
252  * per each particle container and per each centrality bin.
253  */
255 {
256  TString histname;
257  TString histtitle;
258  TString groupname;
259  AliParticleContainer* partCont = 0;
260  TIter next(&fParticleCollArray);
261  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
262  groupname = partCont->GetName();
263  // Protect against creating the histograms twice
264  if (fHistManager.FindObject(groupname)) {
265  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
266  continue;
267  }
268  fHistManager.CreateHistoGroup(groupname);
269  for (Int_t cent = 0; cent < fNcentBins; cent++) {
270  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
271  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
272  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
273 
274  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
275  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
276  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
277 
278  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
279  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
280  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
281 
282  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
283  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
284  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
285  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
286 
287  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
288  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
289  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
290 
291  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
292  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());
293  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
294 
295  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
296  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
297  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
298  }
299 
300  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
301  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
302  if (fForceBeamType != kpp) {
303  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
304  }
305  else {
306  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
307  }
308  }
309  }
310 
311  histname = "fHistSumNTracks";
312  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
313  if (fForceBeamType != kpp) {
314  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
315  }
316  else {
317  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
318  }
319 }
320 
321 /*
322  * This function allocates the histograms for basic jet QA.
323  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
324  * per each jet container and per each centrality bin.
325  */
327 {
328  TString histname;
329  TString histtitle;
330  TString groupname;
331  AliJetContainer* jetCont = 0;
332  TIter next(&fJetCollArray);
333  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
334  groupname = jetCont->GetName();
335  // Protect against creating the histograms twice
336  if (fHistManager.FindObject(groupname)) {
337  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
338  continue;
339  }
340  fHistManager.CreateHistoGroup(groupname);
341  for (Int_t cent = 0; cent < fNcentBins; cent++) {
342  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
343  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
344  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
345 
346  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
347  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
348  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
349 
350  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
351  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
352  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
353 
354  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
355  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
356  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
357 
358  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
359  histtitle = TString::Format("%s;number of jets;events", histname.Data());
360  if (fForceBeamType != kpp) {
361  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
362  }
363  else {
364  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
365  }
366 
367  if (!jetCont->GetRhoName().IsNull()) {
368  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
369  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
370  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
371  }
372  }
373  }
374 }
375 
377 {
378 
379  Bool_t oldStatus = TH1::AddDirectoryStatus();
380  TH1::AddDirectory(kFALSE);
381 
382  // set seed
383  fRandom = new TRandom3(0);
384 
385  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
386  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
387  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
388  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
389  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
390  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
391  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
392 
393  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
394  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
395  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
396  fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
397  fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
398  fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
399  fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
400 
402 
403  fOutput->Add(fh1TrigRef);
404  fOutput->Add(fh1TrigSig);
405  fOutput->Add(fh2Ntriggers);
406  fOutput->Add(fh2RPJetsC10);
407  fOutput->Add(fh2RPJetsC20);
408  fOutput->Add(fh2RPTC10);
409  fOutput->Add(fh2RPTC20);
410 
411  const Int_t dimSpec = 6;
412  const Int_t nBinsSpec[dimSpec] = {100,100, 280, 50,200, fNRPBins};
413  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0,-0.5*TMath::Pi(), 0};
414  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50,1.5*TMath::Pi(), static_cast<Double_t>(fNRPBins)};
415  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
416 
417  // comment out since I want finer binning in jet area, to make it easier
418  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
419  // for inclusive jets)
420  //change binning in jet area
421 // Double_t *xPt6 = new Double_t[7];
422 // xPt6[0] = 0.;
423 // xPt6[1]=0.07;
424 // xPt6[2]=0.2;
425 // xPt6[3]=0.4;
426 // xPt6[4]=0.6;
427 // xPt6[5]=0.8;
428 // xPt6[6]=1;
429 // fHJetSpec->SetBinEdges(1,xPt6);
430 // delete [] xPt6;
431 
432  fOutput->Add(fHJetSpec);
433 
434  // azimuthal correlation
435 
436  fhDphiPtSig = new TH2F("hDphiPtS","recoil #Delta #phi vs jet pT signal",100,-2,5,250,-50,200);
437  fhDphiPtSig->GetXaxis()->SetTitle("#Delta #phi");
438  fhDphiPtSig->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
439  fhDphiPtRef = new TH2F("hDphiPtR","recoil #Delta #phi vs jet pT reference",100,-2,5,250,-50,200);
440  fhDphiPtRef->GetXaxis()->SetTitle("#Delta #phi");
441  fhDphiPtRef->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
442 
443  fOutput->Add(fhDphiPtRef);
444  fOutput->Add(fhDphiPtSig);
445 
446  // =========== Switch on Sumw2 for all histos ===========
447  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
448  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
449  if (h1){
450  h1->Sumw2();
451  continue;
452  }
453  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
454  if (hn){
455  hn->Sumw2();
456  }
457  }
458 
459  // add QA plots from fEventCuts
460  fEventCuts.AddQAplotsToList(fOutput);
461 
462  TH1::AddDirectory(oldStatus);
463 
464  PostData(1, fOutput);
465 }
466 
474 {
475 
476  fHistEvtSelection->Fill(1); // number of events before event selection
477  AliVEvent *ev = InputEvent();
478  if (!fEventCuts.AcceptEvent(ev)) {
479  fHistEvtSelection->Fill(2);
480  return kTRUE;
481  }
482 
485  //DoClusterLoop();
486  //DoCellLoop();
487  DoJetCoreLoop();
488 
489  return kTRUE;
490 }
491 
493 {
494 
495  // Do jet core analysis and fill histograms.
497  if(!jetCont) {
498  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
499  TIter next(&fJetCollArray);
500  while ((jetCont = static_cast<AliJetContainer*>(next())))
501  AliError(Form("%s",jetCont->GetName()));
502  AliFatal("Exit...");
503  return;
504  }
505 
506  // centrality selection
507  if(fDebug) Printf("centrality: %f\n", fCent);
508  if ((fCent>fCentMax) || (fCent<fCentMin)) {
509  fHistEvtSelection->Fill(4);
510  return;
511  }
512  fHistEvtSelection->Fill(0);
513 
514  // Background
515  Double_t rho = 0;
516  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
517  if(fDebug) Printf("rho = %f, rho check = %f",rho, GetRhoVal(0));
518 
519  // Choose trigger track
520  Int_t nT=0;
521  TList ParticleList;
522  Double_t minT=0;
523  Double_t maxT=0;
524  Int_t number=0;
525  Double_t dice=fRandom->Uniform(0,1);
526  Bool_t isSignal = kFALSE;
527  if(dice>fFrac){
528  minT=fTTLowRef;
529  maxT=fTTUpRef;
530  }
531  if(dice<=fFrac){
532  isSignal = kTRUE;
533  minT=fTTLowSig;
534  maxT=fTTUpSig;
535  }
536  nT=SelectTrigger(&ParticleList,minT,maxT,number);
537  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);
538  if(nT<0) return;
539 
540  if(dice>fFrac) fh1TrigRef->Fill(number);
541  if(dice<=fFrac)fh1TrigSig->Fill(number);
542 
543 
544  // particle loop -
545  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
546  // histogram
547  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
548  if(tt!=nT) continue;
549  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
550  if(!partback) continue;
551  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
552  // if(partback->Pt()<8) continue;
553 
554  Int_t injet4=0;
555  Int_t injet=0;
556 
557  fh2Ntriggers->Fill(fCent,partback->Pt());
558  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
559  if(fCent<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
560  if(fCent<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
561 
562  Double_t etabig=0;
563  Double_t ptbig=0;
564  Double_t areabig=0;
565  Double_t phibig=0.;
566  // Double_t areasmall=0;
567 
568  TString histname;
569  TString groupname;
570  groupname = jetCont->GetName();
571  UInt_t count = 0;
572  for(auto jetbig : jetCont->accepted()) {
573  if (!jetbig) continue;
574  count++;
575  ptbig = jetbig->Pt();
576  etabig = jetbig->Eta();
577  phibig = jetbig->Phi();
578  if(ptbig==0) continue;
579  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
580  areabig = jetbig->Area();
581  Double_t ptcorr=ptbig-rho*areabig;
582  //JJJ - perhaps should change eta selection if implemented in jet container
583  if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
584  if(areabig>=0.07) injet=injet+1;
585  if(areabig>=0.4) injet4=injet4+1;
586  Double_t dphi=RelativePhi(partback->Phi(),phibig);
587  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);
588 
589  // do azimuthal correlation analysis
590  // dPhi between -0.5 < dPhi < 1.5
591  Double_t dPhiShift=phibig-partback->Phi();
592  if(dPhiShift>2*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
593  if(dPhiShift<-2*TMath::Pi()) dPhiShift += 2*TMath::Pi();
594  if(dPhiShift<-0.5*TMath::Pi()) dPhiShift += 2*TMath::Pi();
595  if(dPhiShift>1.5*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
596  if(isSignal) fhDphiPtSig->Fill(dPhiShift,ptcorr);
597  else fhDphiPtRef->Fill(dPhiShift,ptcorr);
598 
599  // selection on relative phi
600  if(fJetHadronDeltaPhi>0. &&
601  TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
602 
603  if(fCent<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
604  if(fCent<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
605 
606  Float_t phitt=partback->Phi();
607  if(phitt<0)phitt+=TMath::Pi()*2.;
608  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
609 
610  Double_t fillspec[] = {fCent,areabig,ptcorr,partback->Pt(),dPhiShift, static_cast<Double_t>(phiBintt)};
611  fHJetSpec->Fill(fillspec);
612  }
613  }
614 }
615 
621 {
622  TString histname;
623  TString groupname;
624  AliJetContainer* jetCont = 0;
625  TIter next(&fJetCollArray);
626  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
627  groupname = jetCont->GetName();
628  UInt_t count = 0;
629  for(auto jet : jetCont->accepted()) {
630  if (!jet) continue;
631  count++;
632 
633  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
634  fHistManager.FillTH1(histname, jet->Pt());
635 
636  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
637  fHistManager.FillTH1(histname, jet->Area());
638 
639  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
640  fHistManager.FillTH1(histname, jet->Phi());
641 
642  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
643  fHistManager.FillTH1(histname, jet->Eta());
644 
645  if (jetCont->GetRhoParameter()) {
646  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
647  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
648  }
649  }
650  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
651  fHistManager.FillTH1(histname, count);
652  }
653 }
654 
660 {
662 
663  TString histname;
664  TString groupname;
665  UInt_t sumAcceptedTracks = 0;
666  AliParticleContainer* partCont = 0;
667  TIter next(&fParticleCollArray);
668  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
669  groupname = partCont->GetName();
670  UInt_t count = 0;
671  for(auto part : partCont->accepted()) {
672  if (!part) continue;
673  count++;
674 
675  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
676  fHistManager.FillTH1(histname, part->Pt());
677 
678  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
679  fHistManager.FillTH1(histname, part->Phi());
680 
681  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
682  fHistManager.FillTH1(histname, part->Eta());
683 
684  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
685  const AliVTrack* track = static_cast<const AliVTrack*>(part);
686 
687  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
688  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
689 
690  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
691  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
692 
693  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
694  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
695 
696  if (clusCont) {
697  Int_t iCluster = track->GetEMCALcluster();
698  if (iCluster >= 0) {
699  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
700  if (cluster) {
701  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
702  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
703  }
704  }
705  }
706  }
707  }
708  sumAcceptedTracks += count;
709 
710  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
711  fHistManager.FillTH1(histname, count);
712  }
713 
714  histname = "fHistSumNTracks";
715  fHistManager.FillTH1(histname, sumAcceptedTracks);
716 }
717 
723 {
724  TString histname;
725  TString groupname;
726  UInt_t sumAcceptedClusters = 0;
727  AliClusterContainer* clusCont = 0;
728  TIter next(&fClusterCollArray);
729  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
730  groupname = clusCont->GetName();
731 
732  for(auto cluster : clusCont->all()) {
733  if (!cluster) continue;
734 
735  if (cluster->GetIsExotic()) {
736  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
737  fHistManager.FillTH1(histname, cluster->E());
738  }
739  }
740 
741  UInt_t count = 0;
742  for(auto cluster : clusCont->accepted()) {
743  if (!cluster) continue;
744  count++;
745 
746  AliTLorentzVector nPart;
747  cluster->GetMomentum(nPart, fVertex);
748 
749  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
750  fHistManager.FillTH1(histname, cluster->E());
751 
752  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
753  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
754 
755  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
756  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
757 
758  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
759  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
760 
761  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
762  fHistManager.FillTH1(histname, nPart.Eta());
763  }
764  sumAcceptedClusters += count;
765 
766  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
767  fHistManager.FillTH1(histname, count);
768  }
769 
770  histname = "fHistSumNClusters";
771  fHistManager.FillTH1(histname, sumAcceptedClusters);
772 }
773 
779 {
780  if (!fCaloCells) return;
781 
782  TString histname;
783 
784  const Short_t ncells = fCaloCells->GetNumberOfCells();
785 
786  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
787  fHistManager.FillTH1(histname, ncells);
788 
789  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
790  for (Short_t pos = 0; pos < ncells; pos++) {
791  Double_t amp = fCaloCells->GetAmplitude(pos);
792 
793  fHistManager.FillTH1(histname, amp);
794  }
795 }
796 
802 {
804 }
805 
814 {
815  return kTRUE;
816 }
817 
822 {
823 }
824 
826 
827  Int_t index=-1;
828  Int_t triggers[100];
829 
830  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
831 
832  Int_t im=0;
833 
834  TString groupname = "";
836  groupname = partCont->GetName();
837  UInt_t iCount = 0;
838  for(auto part : partCont->accepted()) {
839  if (!part) continue;
840  list->Add(part);
841  iCount++;
842  if(part->Pt()>=minT && part->Pt()<maxT){
843  triggers[im]=iCount-1;
844  im=im+1;
845  }
846  }
847  number=im;
848  Int_t rd=0;
849  if(im>0) rd=fRandom->Integer(im);
850  index=triggers[rd];
851  return index;
852 }
853 
855 
856  // get relative DeltaPhi from -pi < DeltaPhi < pi
857 
858  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
859  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
860  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
861  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
862  double dphi = mphi-vphi;
863  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
864  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
865  return dphi;//dphi in [-Pi, Pi]
866 }
867 
869 {
870  Int_t phibin=-1;
871  if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
872  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
873  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
874  if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
875  return phibin;
876 }
Float_t fJetHadronDeltaPhi
set 0 > f > pi for selection on jet hadron Dphi
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
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
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 AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Container structure for EMCAL clusters.
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Int_t fNbins
no. of pt bins