AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetMassStructure.cxx
Go to the documentation of this file.
1 //
2 // Jet mass structure analysis task.
3 //
4 // Author: M.Verweij
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <TProfile2D.h>
11 #include <THnSparse.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <TMath.h>
20 #include <TRandom3.h>
21 
22 #include "AliVCluster.h"
23 #include "AliVTrack.h"
24 #include "AliEmcalJet.h"
25 #include "AliRhoParameter.h"
26 #include "AliLog.h"
27 #include "AliEmcalParticle.h"
28 #include "AliMCEvent.h"
29 #include "AliGenPythiaEventHeader.h"
30 #include "AliAODMCHeader.h"
31 #include "AliMCEvent.h"
32 #include "AliAnalysisManager.h"
33 #include "AliParticleContainer.h"
34 #include "AliJetContainer.h"
36 
37 #include "AliAODEvent.h"
38 
40 
42 
43 //________________________________________________________________________
46  fContainerBase(0),
47  fMinFractionShared(0),
48  fJetMassType(kRaw),
49  fRandom(0),
50  fEfficiencyFixed(1.),
51  fCorrType(kMeanPtR),
52  fEJetByJetCorr(0),
53  fh3PtDRMass(0),
54  fh3PtDRRho(0),
55  fh3PtDRMassCorr(0),
56  fh3PtDRRhoCorr(0),
57  fh2PtMass(0),
58  fh2PtMassCorr(0),
59  fhnMassResponse(0),
60  fhnMassResponseCorr(0),
61  fhnDeltaMass(0),
62  fhnDeltaMassCorr(0),
63  fhAllpTRec(0),
64  fhAllpTGen(0),
65  fhAllpTCor(0),
66  fSwitchResolutionhn(1),
67  fh3JetPtDRTrackPt(0),
68  fListOfOutputFromClass(0),
69  fPartArrayN(0)
70 {
71  // Default constructor.
72 
73  fh3PtDRMass = new TH3F*[fNcentBins];
74  fh3PtDRRho = new TH3F*[fNcentBins];
75  fh3PtDRMassCorr = new TH3F*[fNcentBins];
76  fh3PtDRRhoCorr = new TH3F*[fNcentBins];
77  fh2PtMass = new TH2F*[fNcentBins];
78  fh2PtMassCorr = new TH2F*[fNcentBins];
79  fh3JetPtDRTrackPt = new TH3F*[fNcentBins];
80 
81  for (Int_t i = 0; i < fNcentBins; i++) {
82  fh3PtDRMass[i] = 0;
83  fh3PtDRRho[i] = 0;
84  fh3PtDRMassCorr[i] = 0;
85  fh3PtDRRhoCorr[i] = 0;
86  fh2PtMass[i] = 0;
87  fh2PtMassCorr[i] = 0;
88  fh3JetPtDRTrackPt[i] = 0;
89  }
90 
91  SetMakeGeneralHistograms(kTRUE);
92 }
93 
94 //________________________________________________________________________
96  AliAnalysisTaskEmcalJet(name, kTRUE),
97  fContainerBase(0),
98  fMinFractionShared(0),
99  fJetMassType(kRaw),
100  fRandom(0),
101  fEfficiencyFixed(1.),
102  fCorrType(kMeanPtR),
103  fEJetByJetCorr(0),
104  fh3PtDRMass(0),
105  fh3PtDRRho(0),
106  fh3PtDRMassCorr(0),
107  fh3PtDRRhoCorr(0),
108  fh2PtMass(0),
109  fh2PtMassCorr(0),
110  fhnMassResponse(0),
111  fhnMassResponseCorr(0),
112  fhnDeltaMass(0),
113  fhnDeltaMassCorr(0),
114  fhAllpTRec(0),
115  fhAllpTGen(0),
116  fhAllpTCor(0),
117  fSwitchResolutionhn(1),
118  fh3JetPtDRTrackPt(0),
119  fListOfOutputFromClass(0),
120  fPartArrayN(0)
121 {
122  // Standard constructor.
123 
124  fh3PtDRMass = new TH3F*[fNcentBins];
125  fh3PtDRRho = new TH3F*[fNcentBins];
128  fh2PtMass = new TH2F*[fNcentBins];
129  fh2PtMassCorr = new TH2F*[fNcentBins];
131 
132  for (Int_t i = 0; i < fNcentBins; i++) {
133  fh3PtDRMass[i] = 0;
134  fh3PtDRRho[i] = 0;
135  fh3PtDRMassCorr[i] = 0;
136  fh3PtDRRhoCorr[i] = 0;
137  fh2PtMass[i] = 0;
138  fh2PtMassCorr[i] = 0;
139  fh3JetPtDRTrackPt[i] = 0;
140  }
142 }
143 
144 //________________________________________________________________________
146 {
147  // Destructor.
148  delete fRandom;
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create user output.
155 
157 
158  if(!fRandom) fRandom = new TRandom3(0);
159 
161 
162  Bool_t oldStatus = TH1::AddDirectoryStatus();
163  TH1::AddDirectory(kFALSE);
164 
165  const Int_t nBinsPt = 200;
166  const Double_t minPt = -50.;
167  const Double_t maxPt = 150.;
168  Double_t *binsPt = new Double_t[nBinsPt+1];
169  for(Int_t i=0; i<=nBinsPt; i++) binsPt[i]=(Double_t)minPt + (maxPt-minPt)/nBinsPt*(Double_t)i ;
170 
171  const Int_t nBinsM = 200;
172  const Double_t minM = -10.;
173  const Double_t maxM = 40.;
174  Double_t *binsM = new Double_t[nBinsM+1];
175  for(Int_t i=0; i<=nBinsM; i++) binsM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;
176 
177  const Int_t nBinsdM = 200;
178  const Double_t mindM = -30.;
179  const Double_t maxdM = 20.;
180  Double_t *binsdM = new Double_t[nBinsdM+1];
181  for(Int_t i=0; i<=nBinsdM; i++) binsdM[i]=(Double_t)mindM + (maxdM-mindM)/nBinsdM*(Double_t)i ;
182 
183  const Int_t nBinsdMr = 200;
184  const Double_t mindMr = -2.;
185  const Double_t maxdMr = 2.;
186  Double_t *binsdMr = new Double_t[nBinsdMr+1];
187  for(Int_t i=0; i<=nBinsdMr; i++) binsdMr[i]=(Double_t)mindMr + (maxdMr-mindMr)/nBinsdMr*(Double_t)i ;
188 
189  const Int_t nBinsR = 80;
190  const Double_t minR = -0.005;
191  const Double_t maxR = 0.795;
192  Double_t *binsR = new Double_t[nBinsR+1];
193  for(Int_t i=0; i<=nBinsR; i++) binsR[i]=(Double_t)minR + (maxR-minR)/nBinsR*(Double_t)i ;
194 
195  //track pt
196  Double_t binWidth1 = 0.1;
197  Double_t binWidth2 = 1.;
198  Double_t binWidth3 = 1.;
199  const Float_t ptmin1 = 0.;
200  const Float_t ptmax1 = 5.;
201  const Float_t ptmin2 = ptmax1;
202  const Float_t ptmax2 = 10.;
203  const Float_t ptmin3 = ptmax2 ;
204  const Float_t ptmax3 = 100.;
205  const Int_t nbin11 = (int)((ptmax1-ptmin1)/binWidth1);
206  const Int_t nbin12 = (int)((ptmax2-ptmin2)/binWidth2)+nbin11;
207  const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
208  Int_t nBinsPtTr=nbin13;
209  //Create array with low edges of each bin
210  Double_t *binsPtTr=new Double_t[nBinsPtTr+1];
211  for(Int_t i=0; i<=nBinsPtTr; i++) {
212  if(i<=nbin11) binsPtTr[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
213  if(i<=nbin12 && i>nbin11) binsPtTr[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
214  if(i<=nbin13 && i>nbin12) binsPtTr[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
215  }
216 
217  //Binning for THnSparse
218  const Int_t nBinsSparse0 = 4;
219  const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
220  const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt};
221  const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt};
222  const Int_t nBinsSparse1 = 5;
223  const Int_t nBins1[nBinsSparse1] = {nBinsdM,nBinsdMr,nBinsM,nBinsPt,nBinsPt};
224  const Double_t xmin1[nBinsSparse1] = { mindM, mindMr, minM, minPt, minPt};
225  const Double_t xmax1[nBinsSparse1] = { maxdM, maxdMr, maxM, maxPt, maxPt};
226 
227  TString histName = "";
228  TString histTitle = "";
229  for (Int_t i = 0; i < fNcentBins; i++) {
230  histName = TString::Format("fh3PtDRMass_%d",i);
231  histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
232  fh3PtDRMass[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,111,-0.005,1.105);
233  fOutput->Add(fh3PtDRMass[i]);
234 
235  histName = TString::Format("fh3PtDRRho_%d",i);
236  histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
237  fh3PtDRRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,111,-0.005,1.105);
238  fOutput->Add(fh3PtDRRho[i]);
239 
240  histName = TString::Format("fh3PtDRMassCorr_%d",i);
241  histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
242  fh3PtDRMassCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,111,-0.005,1.105);
243  fOutput->Add(fh3PtDRMassCorr[i]);
244 
245  histName = TString::Format("fh3PtDRRhoCorr_%d",i);
246  histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
247  fh3PtDRRhoCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,111,-0.005,1.105);
248  fOutput->Add(fh3PtDRRhoCorr[i]);
249 
250  histName = TString::Format("fh2PtMass_%d",i);
251  histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
252  fh2PtMass[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
253  fOutput->Add(fh2PtMass[i]);
254 
255  histName = TString::Format("fh2PtMassCorr_%d",i);
256  histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
257  fh2PtMassCorr[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
258  fOutput->Add(fh2PtMassCorr[i]);
259 
260  histName = TString::Format("fh3JetPtDRTrackPt_%d",i);
261  histTitle = TString::Format("%s;#it{p}_{T,jet};r;#it{p}_{T,track}",histName.Data());
262  fh3JetPtDRTrackPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,binsPt,nBinsR,binsR,nBinsPtTr,binsPtTr);
263  fOutput->Add(fh3JetPtDRTrackPt[i]);
264  }
265 
266  histName = "fhnMassResponse";
267  histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
268  fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
269  fOutput->Add(fhnMassResponse);
270 
271  histName = "fhnMassResponseCorr";
272  histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
273  fhnMassResponseCorr = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
274  fOutput->Add(fhnMassResponseCorr);
275 
276  if(fEJetByJetCorr) {
277  const Int_t nBinsPtMeanTr = 500;
278  const Double_t minPtMeanTr = 0.;
279  const Double_t maxPtMeanTr = 100.;
280 
285  fhAllpTRec = new TH2F("fhAllpTRec", "Jet constituent p_{T} distribution RECO; #it{p}_{T,track}; #it{p}_{T,jet}", nBinsPtMeanTr, minPtMeanTr, maxPtMeanTr, 20,0.,100.);
286  fhAllpTRec->Sumw2();
287  fhAllpTGen = (TH2F*)fhAllpTRec->Clone("fhAllpTGen");
288  fhAllpTGen->SetTitle("Jet constituent p_{T} distribution GENE; #it{p}_{T,track}; #it{p}_{T,jet}");
289  //fhAllpTCor = new TH2F("fhAllpTCor", "Jet constituent p_{T} distribution CORR; #it{p}_{T,track}; #it{p}_{T,jet}", 500,0.,100., 20,0.,20.);
290  fhAllpTCor = (TH2F*)fhAllpTRec->Clone("fhAllpTCor");
291  fhAllpTCor->SetTitle("Jet constituent p_{T} distribution CORR; #it{p}_{T,track}; #it{p}_{T,jet}");
292 
293  fOutput->Add(fhAllpTRec);
294  fOutput->Add(fhAllpTGen);
295  fOutput->Add(fhAllpTCor);
296 
297  fhtmppTRec = new TH1F("fhtmppTRec", "htmppTRec", 50, 0.,100.);
298  fhtmppTRec->Sumw2();
299  fhtmppTGen = (TH1F*)fhtmppTRec->Clone("fhtmppTGen");
300 
301  const Int_t nvars = 7;
302  Int_t nbins[nvars] = {20 , nBinsPtMeanTr, nBinsPtMeanTr, 20, 20, 50, 50};
303  Double_t limslow[nvars] = {-10, minPtMeanTr, minPtMeanTr, 0., 0., 0., 0.};
304  Double_t limshig[nvars] = {10, maxPtMeanTr, maxPtMeanTr, 100., 100., 1.1, 1.1};
305  TString varnames[nvars] = {"NconstRec-Gen", "pTtmeanRec", "pTtmeanGen", "pTjmeanRec", "pTjmeanGen", "efficiencyNconst", "ConstMatchNormRec"};
306  histName = "fhConstRecGen";
307  histTitle = "Constituents QA; "; for(Int_t i=0;i<nvars; i++) histTitle += Form("%s; ", varnames[i].Data());
308  fhConstRecGen = new THnSparseF(histName, histTitle, nvars, nbins, limslow, limshig);
309  fOutput->Add(fhConstRecGen);
310  }
311 
312  }
314  histName = "fhnDeltaMass";
315  histTitle = Form("%s; #it{M}_{det} - #it{M}_{part} ;(#it{M}_{det} - #it{M}_{part})/#it{M}_{part}; #it{M}_{part}; #it{p}_{T,det}; #it{p}_{T,part}",histName.Data());
316  fhnDeltaMass = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse1,nBins1,xmin1,xmax1);
317  fOutput->Add(fhnDeltaMass);
318 
319  histName = "fhnDeltaMassCorr";
320  histTitle = Form("%s;#it{M}_{det} - #it{M}_{part} ;(#it{M}_{det} - #it{M}_{part})/#it{M}_{part}; #it{M}_{part}; #it{p}_{T,det}; #it{p}_{T,part}",histName.Data());
321  fhnDeltaMassCorr = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse1,nBins1,xmin1,xmax1);
322  fOutput->Add(fhnDeltaMassCorr);
323 
324  }
325 
326  TH1::AddDirectory(oldStatus);
327 
328  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
329 
330  if(binsPt) delete [] binsPt;
331  if(binsPtTr) delete [] binsPtTr;
332  if(binsM) delete [] binsM;
333  if(binsR) delete [] binsR;
334 
335 
336 
337 }
338 
339 //________________________________________________________________________
341 {
342  // Run analysis code here, if needed. It will be executed before FillHistograms().
343 
344  return kTRUE;
345 }
346 
347 //__________________________________________________________________________
349 
351 
352  Int_t Nmissed = 0;
353  TClonesArray* partArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fPartArrayN));
354  if(!partArray){
355  AliError("Array of particles not found, cannot loop into the constituents of the particle level jet");
356  return -1;
357 
358  }
359 
360  if(!jPart || !jet1) return -1;
361 
362  fhtmppTGen->Reset();
363  fhtmppTRec->Reset();
364  Double_t avgpTRec = 0, avgpTGen = 0;
365  Int_t NconstRec = jet1->GetNumberOfTracks(), NconstGen = jPart->GetNumberOfTracks();
366  if(NconstRec == 0) return -1;
367  Int_t labels[NconstRec];
368  for(Int_t i=0 ; i<NconstRec;i++){
369  AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(i, fTracks));
370  if(!vp) continue;
371  Double_t pTtr = vp->Pt();
372  avgpTRec+=pTtr;
373  fhtmppTRec->Fill(pTtr);
374  fhAllpTRec->Fill(pTtr,jet1->Pt());
375  labels[i]=vp->GetLabel();
376  }
377 
378  Int_t countMatch=0;
379  for(Int_t i=0 ; i<NconstGen; i++){
380  AliVParticle *vp = static_cast<AliVParticle*>(jPart->TrackAt(i, partArray));
381  if(!vp) continue;
382  Double_t pTtr = vp->Pt();
383  avgpTGen+=pTtr;
384  fhtmppTGen->Fill(pTtr);
385  fhAllpTGen->Fill(pTtr,jet1->Pt());
386  for(Int_t j=0;j<NconstRec;j++){
387  if (vp->GetLabel()==labels[j]){
388  countMatch++;
389 
390  break;
391  }
392 
393  }
394  }
395 
396  avgpTRec/=(Float_t)NconstRec;
397  avgpTGen/=(Float_t)NconstGen;
398  Double_t fill[7] = {(Double_t)(NconstRec - NconstGen), avgpTRec, avgpTGen, jet1->Pt(), jPart->Pt(), (Float_t)NconstRec/(Float_t)NconstGen, (Double_t)countMatch/(Double_t)NconstRec};
399  fhConstRecGen->Fill(fill);
400  fhtmppTGen->Add(fhtmppTRec, -1);
401  Nmissed = fhtmppTGen->Integral();
402 
403  return Nmissed;
404 
405 }
406 //________________________________________________________________________
408 {
409  // study how jet mass builds up as function of radial distance to jet axis
410 
411  if(fCorrType==kMeanPtR && !fEJetByJetCorr) return kFALSE;
412  AliEmcalJet* jet1 = NULL;
414  Int_t nmissedthisjet = -10;
415  if(jetCont) {
416  jetCont->ResetCurrentID();
417  while((jet1 = jetCont->GetNextAcceptJet())) {
418 
419  Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
420  Double_t mJet1 = GetJetMass(jet1);
421  fh2PtMass[fCentBin]->Fill(ptJet1,mJet1);
422  AliEmcalJet *jPart = jet1->ClosestJet();
423  //fill detector response
424  if(jPart) {
425  nmissedthisjet = CalculateNMissingTracks(jet1, jPart);
426  Double_t var[4] = {mJet1,jPart->M(),ptJet1,jPart->Pt()};
427  fhnMassResponse->Fill(var);
428  Double_t var1[5] = {mJet1-jPart->M(),(mJet1-jPart->M())/jPart->M(),jPart->M(), ptJet1,jPart->Pt()};
429  fhnDeltaMass->Fill(var1);
430 
431  }
432 
433  // if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
434  // continue;
435 
436  //Sort array based on distance to jet axis
437  static Int_t indexes[9999] = {-1};
438  static Float_t dr[9999] = {0};
439  if(!GetSortedArray(indexes,dr,jet1)) continue;
440 
441  for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
442  AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
443  if(!vp) continue;
444  fh3JetPtDRTrackPt[fCentBin]->Fill(ptJet1,dr[indexes[i]],vp->Pt());
445  }
446 
447  if(fCorrType==kAnnulus) {
448  TLorentzVector sumVec; sumVec.SetPtEtaPhiM(0.,0.,0.,0.);
449  TLorentzVector corrVec; corrVec.SetPtEtaPhiM(0.,0.,0.,0.);
450  TLorentzVector curVec;
451  for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
452  AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
453  if(!vp) continue;
454 
455  curVec.SetPtEtaPhiM(vp->Pt(),vp->Eta(),vp->Phi(),vp->M());
456  sumVec+=curVec;
457  corrVec+=curVec;
458 
459  fh3PtDRMass[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.M()/mJet1);
460  fh3PtDRRho[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.Pt()/ptJet1);
461 
462  fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
463  fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
464 
465  Double_t eff = GetEfficiency(vp->Pt());
466  Double_t rnd = fRandom->Uniform(1.);
467  if(rnd>eff) {//put particle back at random position in annulus; using now zero width for annulus
468  Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
469  Double_t rr = dr[indexes[i]];
470  rr = fRandom->Uniform(0.,jetCont->GetJetRadius());
471  Double_t deta = rr*TMath::Cos(t);
472  Double_t dphi = rr*TMath::Sin(t);
473  curVec.SetPtEtaPhiM(vp->Pt(),deta+jet1->Eta(),dphi+jet1->Phi(),vp->M());
474  corrVec+=curVec;
475 
476  fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
477  fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
478  }
479  }//track loop
480 
481  fh2PtMassCorr[fCentBin]->Fill(corrVec.Pt(), corrVec.M());
482  if(jPart) {
483  Double_t varCorr[4] = {corrVec.M(),jPart->M(),corrVec.Pt(),jPart->Pt()};
484  fhnMassResponseCorr->Fill(varCorr);
485  }
486  }//kAnnulus method
487  else if(fCorrType==kMeanPtR) {
488  //jet-by-jet correction based on templates
489  if(nmissedthisjet > 0)
490  fEJetByJetCorr->SetNMissedTracks(nmissedthisjet);
491 
492  AliEmcalJet *jetCorr = fEJetByJetCorr->Eval(jet1,jetCont->GetParticleContainer()->GetArray());
493  if(jPart && jetCorr) {
494  Double_t varCorr[4] = {jetCorr->M(),jPart->M(),jetCorr->Pt(),jPart->Pt()};
495  fhnMassResponseCorr->Fill(varCorr);
496  Double_t varCorr1[5] = {jetCorr->M()-jPart->M(),(jetCorr->M()-jPart->M())/jPart->M(),jPart->M(),jetCorr->Pt(),jPart->Pt()};
497  fhnDeltaMassCorr->Fill(varCorr1);
499 
500  TClonesArray* partArray = fEJetByJetCorr->GetAddedTrackArray();
501  if(!partArray){
502  AliError("Array of added particles not found, cannot loop into the constituents");
503  return kFALSE;
504 
505  }
506 
507  Int_t NconstCorr = partArray->GetEntries();
508 
509  for(Int_t l=0 ; l<NconstCorr; l++){
510  TLorentzVector *vp = (TLorentzVector*)partArray->At(l);
511  if(!vp) continue;
512  Double_t pTtr = vp->Pt();
513  fhAllpTCor->Fill(pTtr,jet1->Pt());
514 
515  }
516 
517  }
518  }
519  }//kMeanPtR method
520  nmissedthisjet = -10;
521 
522  }//jet loop
523  }//jet container
524 
525 
526 
527  return kTRUE;
528 }
529 
530 //________________________________________________________________________
532  pt = 1.*pt;
533  Double_t eff = 1.;
534  if(fEfficiencyFixed<1.) return fEfficiencyFixed;
535 
536  return eff;
537 }
538 
539 //________________________________________________________________________
541 {
542  // sort array
543  const Int_t n = (Int_t)jet->GetNumberOfTracks();//(Int_t)array.size();
544  if (n < 1) return kFALSE;
545 
546  for (Int_t i = 0; i < n; i++) {
547  AliVParticle *vp = static_cast<AliVParticle*>(jet->TrackAt(i, fTracks));
548  if(!vp) continue;
549  dr[i] = jet->DeltaR(vp);
550  }
551  TMath::Sort(n, dr, indexes);
552  return kTRUE;
553 }
554 
555 
556 //________________________________________________________________________
558  //calc subtracted jet mass
559  if(fJetMassType==kRaw)
560  return jet->M();
561  else if(fJetMassType==kDeriv)
563 
564  return 0;
565 }
566 
567 //________________________________________________________________________
569  //
570  // retrieve event objects
571  //
573  return kFALSE;
574 
575  return kTRUE;
576 }
577 
578 //_______________________________________________________________________
580 {
581  // Called once at the end of the analysis.
582 }
583 
TString fPartArrayN
Array of particles used for jet reconstruction at particle level (need to make it transient probably)...
Double_t Area() const
Definition: AliEmcalJet.h:123
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:219
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Phi() const
Definition: AliEmcalJet.h:110
TH3F ** fh3PtDRMassCorr
! jet pT vs dr(jet axis, constituent) vs cumulative mass density corrected
THnSparse * fhnMassResponseCorr
! response matrix corrected
TList * fListOfOutputFromClass
! list of output from class AliEmcalJetByJetCorrection
AliEmcalJet * Eval(const AliEmcalJet *jet, TClonesArray *fTracks)
void SetNMissedTracks(Double_t number)
THnSparse * fhnDeltaMass
! resolution on mass matrix
TH2F * fhAllpTRec
! histogram that stores the pT of the constituents (RECO level)
TH2F * fhAllpTGen
! histogram that stores the pT of the constituents (PART level)
Int_t fCentBin
!event centrality bin
TRandom * gRandom
THnSparse * fhnDeltaMassCorr
! resolution on mass matrix corrected
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
TH3F ** fh3JetPtDRTrackPt
! jet pt vs dr(jet axis, constituent) vs pT,track
AliParticleContainer * GetParticleContainer() const
JetByJetCorrType fCorrType
jet-by-jet correction method
Double_t fEfficiencyFixed
fixed efficiency for all pT and all types of tracks
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
TH3F ** fh3PtDRRho
! jet pT vs dr(jet axis, constituent) vs cumulative pt density
Int_t fNcentBins
how many centrality bins
TH3F ** fh3PtDRRhoCorr
! jet pT vs dr(jet axis, constituent) vs cumulative pt density corrected
Bool_t GetSortedArray(Int_t indexes[], Float_t dr[], AliEmcalJet *jet) const
AliEmcalJet * GetNextAcceptJet()
Double_t DeltaR(const AliVParticle *part) const
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t Pt() const
Definition: AliEmcalJet.h:102
JetMassType fJetMassType
jet mass type to be used
Double_t GetRhoVal(Int_t i=0) const
TH3F ** fh3PtDRMass
! jet pT vs dr(jet axis, constituent) vs cumulative mass density
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
TClonesArray * fTracks
!tracks
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
AliEmcalJetByJetCorrection * fEJetByJetCorr
object to do jet-by-jet correction
Bool_t fSwitchResolutionhn
switch on/off (default on) the 2 THnSparse for the mass resolution
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48
const Int_t nbins
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:253
TH2F * fhAllpTCor
! histogram that stores the pT of the constituents (CORR level)
Int_t CalculateNMissingTracks(AliEmcalJet *jet1, AliEmcalJet *jPart)
TClonesArray * GetAddedTrackArray() const
Double_t M() const
Definition: AliEmcalJet.h:113
ClassImp(AliAnalysisTaskEmcalJetMassStructure) AliAnalysisTaskEmcalJetMassStructure
TH1F * fhtmppTGen
Temporary stores the pT distribution of jet constituents (gen level)
Container for jet within the EMCAL jet framework.
THnSparse * fhConstRecGen
! number of constituent correlation
TH1F * fhtmppTRec
Temporary stores the pT distribution of jet constituents (reco level)