AliPhysics  8d00e07 (8d00e07)
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  fRunAnaAzimuthalCorrelation(kFALSE),
68  fRandom(0),
69  fHistEvtSelection(0x0),
70  fHJetSpec(0x0),
71  fh1TrigRef(0x0),
72  fh1TrigSig(0x0),
73  fh2Ntriggers(0x0),
74  fh2RPJetsC10(0x0),
75  fh2RPJetsC20(0x0),
76  fh2RPTC10(0x0),
77  fh2RPTC20(0x0),
78  fHJetPhiCorr(0x0),
79  fhDphiPtSig(0x0),
80  fhDphiPtRef(0x0)
81 {
82 }
83 
90  AliAnalysisTaskEmcalJet(name, kTRUE),
91  fHistManager(name),
92  fCentMin(0.),
93  fCentMax(100.),
94  fTTLowRef(8),
95  fTTUpRef(9.),
96  fTTLowSig(20.),
97  fTTUpSig(50.),
98  fNRPBins(50),
99  fFrac(0.8),
100  fJetEtaMin(-.5),
101  fJetEtaMax(.5),
102  fJetHadronDeltaPhi(0.6),
103  fJetContName(""),
105  fRandom(0),
106  fHistEvtSelection(0x0),
107  fHJetSpec(0x0),
108  fh1TrigRef(0x0),
109  fh1TrigSig(0x0),
110  fh2Ntriggers(0x0),
111  fh2RPJetsC10(0x0),
112  fh2RPJetsC20(0x0),
113  fh2RPTC10(0x0),
114  fh2RPTC20(0x0),
115  fHJetPhiCorr(0x0),
116  fhDphiPtSig(0x0),
117  fhDphiPtRef(0x0)
118 {
120 }
121 
126 {
127 }
128 
134 {
136 
142 
143  TIter next(fHistManager.GetListOfHistograms());
144  TObject* obj = 0;
145  while ((obj = next())) {
146  fOutput->Add(obj);
147  }
148 
149  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
150 }
151 
152 /*
153  * This function allocates the histograms for basic EMCal cluster QA.
154  * A set of histograms (energy, eta, phi, number of cluster) is allocated
155  * per each cluster container and per each centrality bin.
156  */
158 {
159  TString histname;
160  TString histtitle;
161  TString groupname;
162  AliClusterContainer* clusCont = 0;
163  TIter next(&fClusterCollArray);
164  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
165  groupname = clusCont->GetName();
166  // Protect against creating the histograms twice
167  if (fHistManager.FindObject(groupname)) {
168  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
169  continue;
170  }
171  fHistManager.CreateHistoGroup(groupname);
172  for (Int_t cent = 0; cent < fNcentBins; cent++) {
173  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
174  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
175  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
176 
177  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
178  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
179  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
180 
181  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
182  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
183  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
184 
185  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
186  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
187  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
188 
189  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
190  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
191  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
192 
193  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
194  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
195  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
196 
197  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
198  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
199  if (fForceBeamType != kpp) {
200  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
201  }
202  else {
203  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
204  }
205  }
206  }
207 
208  histname = "fHistSumNClusters";
209  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
210  if (fForceBeamType != kpp) {
211  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
212  }
213  else {
214  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
215  }
216 }
217 
218 /*
219  * This function allocates the histograms for basic EMCal QA.
220  * One 2D histogram with the cell energy spectra and the number of cells
221  * per event is allocated per each centrality bin.
222  */
224 {
225  TString histname;
226  TString histtitle;
227  TString groupname(fCaloCellsName);
228 
229  fHistManager.CreateHistoGroup(groupname);
230  for (Int_t cent = 0; cent < fNcentBins; cent++) {
231  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
232  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
233  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
234 
235  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
236  histtitle = TString::Format("%s;number of cells;events", histname.Data());
237  if (fForceBeamType != kpp) {
238  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
239  }
240  else {
241  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
242  }
243  }
244 }
245 
246 /*
247  * This function allocates the histograms for basic tracking QA.
248  * A set of histograms (pT, eta, phi, difference between kinematic properties
249  * at the vertex and at the EMCal surface, number of tracks) is allocated
250  * per each particle container and per each centrality bin.
251  */
253 {
254  TString histname;
255  TString histtitle;
256  TString groupname;
257  AliParticleContainer* partCont = 0;
258  TIter next(&fParticleCollArray);
259  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
260  groupname = partCont->GetName();
261  // Protect against creating the histograms twice
262  if (fHistManager.FindObject(groupname)) {
263  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
264  continue;
265  }
266  fHistManager.CreateHistoGroup(groupname);
267  for (Int_t cent = 0; cent < fNcentBins; cent++) {
268  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
269  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
270  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
271 
272  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
273  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
274  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
275 
276  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
277  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
278  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
279 
280  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
281  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
282  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
283  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
284 
285  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
286  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
287  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
288 
289  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
290  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());
291  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
292 
293  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
294  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
295  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
296  }
297 
298  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
299  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
300  if (fForceBeamType != kpp) {
301  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
302  }
303  else {
304  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
305  }
306  }
307  }
308 
309  histname = "fHistSumNTracks";
310  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
311  if (fForceBeamType != kpp) {
312  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
313  }
314  else {
315  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
316  }
317 }
318 
319 /*
320  * This function allocates the histograms for basic jet QA.
321  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
322  * per each jet container and per each centrality bin.
323  */
325 {
326  TString histname;
327  TString histtitle;
328  TString groupname;
329  AliJetContainer* jetCont = 0;
330  TIter next(&fJetCollArray);
331  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
332  groupname = jetCont->GetName();
333  // Protect against creating the histograms twice
334  if (fHistManager.FindObject(groupname)) {
335  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
336  continue;
337  }
338  fHistManager.CreateHistoGroup(groupname);
339  for (Int_t cent = 0; cent < fNcentBins; cent++) {
340  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
341  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
342  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
343 
344  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
345  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
346  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
347 
348  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
349  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
350  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
351 
352  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
353  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
354  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
355 
356  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
357  histtitle = TString::Format("%s;number of jets;events", histname.Data());
358  if (fForceBeamType != kpp) {
359  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
360  }
361  else {
362  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
363  }
364 
365  if (!jetCont->GetRhoName().IsNull()) {
366  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
367  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
368  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
369  }
370  }
371  }
372 }
373 
375 {
376 
377  Bool_t oldStatus = TH1::AddDirectoryStatus();
378  TH1::AddDirectory(kFALSE);
379 
380  // set seed
381  fRandom = new TRandom3(0);
382 
383  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
384  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
385  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
386  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
387  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
388  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
389  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
390 
391  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
392  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
393  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
394  fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
395  fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
396  fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
397  fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
398 
400 
401  fOutput->Add(fh1TrigRef);
402  fOutput->Add(fh1TrigSig);
403  fOutput->Add(fh2Ntriggers);
404  fOutput->Add(fh2RPJetsC10);
405  fOutput->Add(fh2RPJetsC20);
406  fOutput->Add(fh2RPTC10);
407  fOutput->Add(fh2RPTC20);
408 
409  const Int_t dimSpec = 5;
410  const Int_t nBinsSpec[dimSpec] = {100,100, 140, 50, fNRPBins};
411  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0, 0};
412  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50, static_cast<Double_t>(fNRPBins)};
413  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
414 
415  // comment out since I want finer binning in jet area, to make it easier
416  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
417  // for inclusive jets)
418  //change binning in jet area
419 // Double_t *xPt6 = new Double_t[7];
420 // xPt6[0] = 0.;
421 // xPt6[1]=0.07;
422 // xPt6[2]=0.2;
423 // xPt6[3]=0.4;
424 // xPt6[4]=0.6;
425 // xPt6[5]=0.8;
426 // xPt6[6]=1;
427 // fHJetSpec->SetBinEdges(1,xPt6);
428 // delete [] xPt6;
429 
430  fOutput->Add(fHJetSpec);
431 
432  // azimuthal correlation
434 
435  const Int_t dimCor = 5;
436  const Int_t nBinsCor[dimCor] = {50, 200, 100, 100, 100};
437  const Double_t lowBinCor[dimCor] = {0, -50, -0.5*TMath::Pi(), 0, 0};
438  const Double_t hiBinCor[dimCor] = {50, 150, 1.5*TMath::Pi(), 1, 100};
439  fHJetPhiCorr = new THnSparseF("fHJetPhiCorr","TT p_{T} vs jet p_{T} vs dPhi vs area vs centrality",dimCor,nBinsCor,lowBinCor,hiBinCor);
440 
441  fhDphiPtSig = new TH2F("hDphiPtS","recoil #Delta #phi vs jet pT signal",100,-2,5,250,-50,200);
442  fhDphiPtSig->GetXaxis()->SetTitle("#Delta #phi");
443  fhDphiPtSig->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
444  fhDphiPtRef = new TH2F("hDphiPtR","recoil #Delta #phi vs jet pT reference",100,-2,5,250,-50,200);
445  fhDphiPtRef->GetXaxis()->SetTitle("#Delta #phi");
446  fhDphiPtRef->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
447 
448  fOutput->Add(fHJetPhiCorr);
449  fOutput->Add(fhDphiPtRef);
450  fOutput->Add(fhDphiPtSig);
451  }
452 
453  // =========== Switch on Sumw2 for all histos ===========
454  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
455  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
456  if (h1){
457  h1->Sumw2();
458  continue;
459  }
460  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
461  if (hn){
462  hn->Sumw2();
463  }
464  }
465 
466  // add QA plots from fEventCuts
467  fEventCuts.AddQAplotsToList(fOutput);
468 
469  TH1::AddDirectory(oldStatus);
470 
471  PostData(1, fOutput);
472 }
473 
481 {
482 
483  fHistEvtSelection->Fill(1); // number of events before event selection
484  AliVEvent *ev = InputEvent();
485  if (!fEventCuts.AcceptEvent(ev)) {
486  fHistEvtSelection->Fill(2);
487  return kTRUE;
488  }
489 
490  DoJetLoop();
491  DoTrackLoop();
492  //DoClusterLoop();
493  //DoCellLoop();
494  DoJetCoreLoop();
495 
496  return kTRUE;
497 }
498 
500 {
501 
502  // Do jet core analysis and fill histograms.
504  if(!jetCont) {
505  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
506  TIter next(&fJetCollArray);
507  while ((jetCont = static_cast<AliJetContainer*>(next())))
508  AliError(Form("%s",jetCont->GetName()));
509  AliFatal("Exit...");
510  return;
511  }
512 
513  // centrality selection
514  if(fDebug) Printf("centrality: %f\n", fCent);
515  if ((fCent>fCentMax) || (fCent<fCentMin)) {
516  fHistEvtSelection->Fill(4);
517  return;
518  }
519  fHistEvtSelection->Fill(0);
520 
521  // Background
522  Double_t rho = 0;
523  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
524  if(fDebug) Printf("rho = %f, rho check = %f",rho, GetRhoVal(0));
525 
526  // Choose trigger track
527  Int_t nT=0;
528  TList ParticleList;
529  Double_t minT=0;
530  Double_t maxT=0;
531  Int_t number=0;
532  Double_t dice=fRandom->Uniform(0,1);
533  Bool_t isSignal = kFALSE;
534  if(dice>fFrac){
535  minT=fTTLowRef;
536  maxT=fTTUpRef;
537  }
538  if(dice<=fFrac){
539  isSignal = kTRUE;
540  minT=fTTLowSig;
541  maxT=fTTUpSig;
542  }
543  nT=SelectTrigger(&ParticleList,minT,maxT,number);
544  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);
545  if(nT<0) return;
546 
547  if(dice>fFrac) fh1TrigRef->Fill(number);
548  if(dice<=fFrac)fh1TrigSig->Fill(number);
549 
550 
551  // particle loop -
552  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
553  // histogram
554  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
555  if(tt!=nT) continue;
556  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
557  if(!partback) continue;
558  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
559  // if(partback->Pt()<8) continue;
560 
561  Int_t injet4=0;
562  Int_t injet=0;
563 
564  fh2Ntriggers->Fill(fCent,partback->Pt());
565  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
566  if(fCent<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
567  if(fCent<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
568 
569  Double_t etabig=0;
570  Double_t ptbig=0;
571  Double_t areabig=0;
572  Double_t phibig=0.;
573  // Double_t areasmall=0;
574 
575  TString histname;
576  TString groupname;
577  groupname = jetCont->GetName();
578  UInt_t count = 0;
579  for(auto jetbig : jetCont->accepted()) {
580  if (!jetbig) continue;
581  count++;
582  ptbig = jetbig->Pt();
583  etabig = jetbig->Eta();
584  phibig = jetbig->Phi();
585  if(ptbig==0) continue;
586  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
587  areabig = jetbig->Area();
588  Double_t ptcorr=ptbig-rho*areabig;
589  //JJJ - perhaps should change eta selection if implemented in jet container
590  if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
591  if(areabig>=0.07) injet=injet+1;
592  if(areabig>=0.4) injet4=injet4+1;
593  Double_t dphi=RelativePhi(partback->Phi(),phibig);
594  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);
595 
596  // do azimuthal correlation analysis
598 
599  // dPhi between -0.5 < dPhi < 1.5
600  Double_t dPhiShift=phibig-partback->Phi();
601  if(dPhiShift>2*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
602  if(dPhiShift<-2*TMath::Pi()) dPhiShift += 2*TMath::Pi();
603  if(dPhiShift<-0.5*TMath::Pi()) dPhiShift += 2*TMath::Pi();
604  if(dPhiShift>1.5*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
605  if(isSignal) fhDphiPtSig->Fill(dPhiShift,ptcorr);
606  else fhDphiPtRef->Fill(dPhiShift,ptcorr);
607  Double_t fill[] = {partback->Pt(),ptcorr,dPhiShift,areabig,fCent};
608  fHJetPhiCorr->Fill(fill);
609  }
610  // selection on relative phi
611  if(TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
612 
613  if(fCent<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
614  if(fCent<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
615 
616  Float_t phitt=partback->Phi();
617  if(phitt<0)phitt+=TMath::Pi()*2.;
618  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
619 
620  Double_t fillspec[] = {fCent,jetbig->Area(),ptcorr,partback->Pt(), static_cast<Double_t>(phiBintt)};
621  fHJetSpec->Fill(fillspec);
622  }
623 
624  //Implementation in old task
625 // if(fRunAnaAzimuthalCorrelation) {
626 // for(auto jetbig : jetCont->accepted()) {
627 // if (!jetbig) continue;
628 // Double_t jetPt = jetbig->Pt();
629 // Double_t jetEta = jetbig->Eta();
630 // Double_t jetPhi = jetbig->Phi();
631 // if(jetPt==0) continue;
632 // if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
633 // Double_t jetArea = jetbig->EffectiveAreaCharged();
634 // Double_t jetPtCorr=jetPt-rho*jetArea;
635 // Double_t dPhi=jetPhi-partback->Phi();
636 // if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
637 // if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
638 // if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
639 // if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
640 // if(fDebug) Printf("\t phi az = %f \tphi function = %f",dPhi,dphi);
641 // Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
642 // fHJetPhiCorr->Fill(fill);
643 // }
644 // }
645 
646 
647 
648  }
649 }
650 
656 {
657  TString histname;
658  TString groupname;
659  AliJetContainer* jetCont = 0;
660  TIter next(&fJetCollArray);
661  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
662  groupname = jetCont->GetName();
663  UInt_t count = 0;
664  for(auto jet : jetCont->accepted()) {
665  if (!jet) continue;
666  count++;
667 
668  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
669  fHistManager.FillTH1(histname, jet->Pt());
670 
671  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
672  fHistManager.FillTH1(histname, jet->Area());
673 
674  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
675  fHistManager.FillTH1(histname, jet->Phi());
676 
677  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
678  fHistManager.FillTH1(histname, jet->Eta());
679 
680  if (jetCont->GetRhoParameter()) {
681  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
682  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
683  }
684  }
685  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
686  fHistManager.FillTH1(histname, count);
687  }
688 }
689 
695 {
697 
698  TString histname;
699  TString groupname;
700  UInt_t sumAcceptedTracks = 0;
701  AliParticleContainer* partCont = 0;
702  TIter next(&fParticleCollArray);
703  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
704  groupname = partCont->GetName();
705  UInt_t count = 0;
706  for(auto part : partCont->accepted()) {
707  if (!part) continue;
708  count++;
709 
710  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
711  fHistManager.FillTH1(histname, part->Pt());
712 
713  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
714  fHistManager.FillTH1(histname, part->Phi());
715 
716  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
717  fHistManager.FillTH1(histname, part->Eta());
718 
719  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
720  const AliVTrack* track = static_cast<const AliVTrack*>(part);
721 
722  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
723  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
724 
725  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
726  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
727 
728  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
729  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
730 
731  if (clusCont) {
732  Int_t iCluster = track->GetEMCALcluster();
733  if (iCluster >= 0) {
734  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
735  if (cluster) {
736  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
737  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
738  }
739  }
740  }
741  }
742  }
743  sumAcceptedTracks += count;
744 
745  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
746  fHistManager.FillTH1(histname, count);
747  }
748 
749  histname = "fHistSumNTracks";
750  fHistManager.FillTH1(histname, sumAcceptedTracks);
751 }
752 
758 {
759  TString histname;
760  TString groupname;
761  UInt_t sumAcceptedClusters = 0;
762  AliClusterContainer* clusCont = 0;
763  TIter next(&fClusterCollArray);
764  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
765  groupname = clusCont->GetName();
766 
767  for(auto cluster : clusCont->all()) {
768  if (!cluster) continue;
769 
770  if (cluster->GetIsExotic()) {
771  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
772  fHistManager.FillTH1(histname, cluster->E());
773  }
774  }
775 
776  UInt_t count = 0;
777  for(auto cluster : clusCont->accepted()) {
778  if (!cluster) continue;
779  count++;
780 
781  AliTLorentzVector nPart;
782  cluster->GetMomentum(nPart, fVertex);
783 
784  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
785  fHistManager.FillTH1(histname, cluster->E());
786 
787  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
788  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
789 
790  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
791  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
792 
793  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
794  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
795 
796  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
797  fHistManager.FillTH1(histname, nPart.Eta());
798  }
799  sumAcceptedClusters += count;
800 
801  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
802  fHistManager.FillTH1(histname, count);
803  }
804 
805  histname = "fHistSumNClusters";
806  fHistManager.FillTH1(histname, sumAcceptedClusters);
807 }
808 
814 {
815  if (!fCaloCells) return;
816 
817  TString histname;
818 
819  const Short_t ncells = fCaloCells->GetNumberOfCells();
820 
821  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
822  fHistManager.FillTH1(histname, ncells);
823 
824  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
825  for (Short_t pos = 0; pos < ncells; pos++) {
826  Double_t amp = fCaloCells->GetAmplitude(pos);
827 
828  fHistManager.FillTH1(histname, amp);
829  }
830 }
831 
837 {
839 }
840 
849 {
850  return kTRUE;
851 }
852 
857 {
858 }
859 
860 THnSparse* AliAnalysisTaskJetCoreEmcal::NewTHnSparseF(const char* name, UInt_t entries)
861 {
862  // generate new THnSparseF, axes are defined in GetDimParams()
863 
864  Int_t count = 0;
865  UInt_t tmp = entries;
866  while(tmp!=0){
867  count++;
868  tmp = tmp &~ -tmp; // clear lowest bit
869  }
870 
871  TString hnTitle(name);
872  const Int_t dim = count;
873  Int_t nbins[dim];
874  Double_t xmin[dim];
875  Double_t xmax[dim];
876 
877  Int_t i=0;
878  Int_t c=0;
879  while(c<dim && i<32){
880  if(entries&(1<<i)){
881 
882  TString label("");
883  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
884  hnTitle += Form(";%s",label.Data());
885  c++;
886  }
887 
888  i++;
889  }
890  hnTitle += ";";
891 
892  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
893 }
894 
896 {
897 
898  // stores label and binning of axis for THnSparse
899 
900  const Double_t pi = TMath::Pi();
901  switch(iEntry){
902  case 0:
903  label = "V0 centrality (%)";
904  nbins = 10;
905  xmin = 0.;
906  xmax = 100.;
907  break;
908  case 1:
909  label = "corrected jet pt";
910  nbins = 20;
911  xmin = 0.;
912  xmax = 200.;
913  break;
914  case 2:
915  label = "track pT";
916  nbins = 9;
917  xmin = 0.;
918  xmax = 150;
919  break;
920  case 3:
921  label = "deltaR";
922  nbins = 15;
923  xmin = 0.;
924  xmax = 1.5;
925  break;
926  case 4:
927  label = "deltaEta";
928  nbins = 8;
929  xmin = -1.6;
930  xmax = 1.6;
931  break;
932  case 5:
933  label = "deltaPhi";
934  nbins = 90;
935  xmin = -0.5*pi;
936  xmax = 1.5*pi;
937  break;
938  case 6:
939  label = "leading track";
940  nbins = 13;
941  xmin = 0;
942  xmax = 50;
943  break;
944  case 7:
945  label = "trigger track";
946  nbins =10;
947  xmin = 0;
948  xmax = 50;
949  break;
950  }
951 }
952 
954 
955  Int_t index=-1;
956  Int_t triggers[100];
957 
958  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
959 
960  Int_t im=0;
961 
962  TString groupname = "";
964  groupname = partCont->GetName();
965  UInt_t iCount = 0;
966  for(auto part : partCont->accepted()) {
967  if (!part) continue;
968  list->Add(part);
969  iCount++;
970  if(part->Pt()>=minT && part->Pt()<maxT){
971  triggers[im]=iCount-1;
972  im=im+1;
973  }
974  }
975  number=im;
976  Int_t rd=0;
977  if(im>0) rd=fRandom->Integer(im);
978  index=triggers[rd];
979  return index;
980 }
981 
983 
984  // get relative DeltaPhi from -pi < DeltaPhi < pi
985 
986  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
987  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
988  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
989  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
990  double dphi = mphi-vphi;
991  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
992  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
993  return dphi;//dphi in [-Pi, Pi]
994 }
995 
997 {
998  Int_t phibin=-1;
999  if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1000  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1001  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1002  if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1003  return phibin;
1004 }
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
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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