AliPhysics  6bc8652 (6bc8652)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalJetByJetCorrection.cxx
Go to the documentation of this file.
1 //
2 // Utility class to do jet-by-jet correction
3 // Based on templates of ptJet vs ptTrack vs r
4 //
5 // Author: M.Verweij
6 
7 #include "TH2.h"
8 #include "TH3.h"
9 #include "TProfile.h"
10 #include "TMath.h"
11 #include "TRandom3.h"
12 #include "TLorentzVector.h"
13 #include "AliEmcalJet.h"
14 #include "AliVParticle.h"
15 #include <TF1.h>
16 #include <TList.h>
17 #include <THnSparse.h>
18 
20 
22 
23 //__________________________________________________________________________
25 TNamed(),
26  fh3JetPtDRTrackPt(0x0),
27  fBinWidthJetPt(10.),
28  fJetPtMin(0.),
29  fJetPtMax(150.),
30  fCollTemplates(),
31  fInitialized(kFALSE),
32  fEfficiencyFixed(1.),
33  fhEfficiency(0),
34  fhSmoothEfficiency(0),
35  fCorrectpTtrack(0),
36  fNpPoisson(0),
37  fExternalNmissed(0),
38  fRndm(0),
39  fNMissedTracks(-1),
40  fpAppliedEfficiency(0),
41  fhNmissing(0),
42  fListOfOutput(0)
43 {
44  // Dummy constructor.
45  fCollTemplates.SetOwner(kTRUE);
46 }
47 
48 //__________________________________________________________________________
50  TNamed(name, name),
51  fh3JetPtDRTrackPt(0x0),
52  fBinWidthJetPt(10.),
53  fJetPtMin(0.),
54  fJetPtMax(150.),
55  fCollTemplates(),
56  fInitialized(kFALSE),
57  fEfficiencyFixed(1.),
58  fhEfficiency(0),
59  fhSmoothEfficiency(0),
60  fCorrectpTtrack(0),
61  fNpPoisson(0),
62  fExternalNmissed(0),
63  fNMissedTracks(-1),
64  fRndm(0),
65  fpAppliedEfficiency(0),
66  fhNmissing(0),
67  fListOfOutput(0)
68 {
69  // Default constructor.
70  fCollTemplates.SetOwner(kTRUE);
71 
72  const Int_t nBinPt = 118; //0-2: 20 bins; 2-100: 98 bins
73  Double_t binLimitsPt[nBinPt+1];
74  for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
75  if(iPt<20){
76  binLimitsPt[iPt] = 0. + (Double_t)iPt*0.15;
77  } else {// 1.0
78  binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
79  }
80  }
81  const Int_t nBinsPtJ = 200;
82  const Double_t minPtJ = -50.;
83  const Double_t maxPtJ = 150.;
84 
85  fpAppliedEfficiency = new TProfile("fpAppliedEfficiency","fpAppliedEfficiency",nBinPt,binLimitsPt);
86 
87  const Int_t nvars = 5;
88  Int_t nbins[nvars] = {nBinsPtJ, 21 , 21 , 21 , 21};
89  Double_t minbin[nvars] = {minPtJ , 0. , 0. , 0. , 0.};
90  Double_t maxbin[nvars] = {maxPtJ , 20., 20., 20., 20.};
91  TString title = "fhNmissing", nameh = title;
92  TString axtitles[nvars] = {"#it{p}_{T,jet}", "N constituents added", "N_{constituents} #times (1/eff - 1)", "N_{truth}"};
93  for(Int_t i = 0; i<nvars; i++){
94  title+=axtitles[i];
95  }
96  fhNmissing = new THnSparseF(nameh.Data(), title.Data(), nvars, nbins, minbin, maxbin);
97 
98  fListOfOutput = new TList();
99  fListOfOutput->SetName("JetByJetCorrectionOutput");
100  fListOfOutput->SetOwner();
101  fListOfOutput->Add(fpAppliedEfficiency);
102  fListOfOutput->Add(fhNmissing);
103 
104 
105 }
106 
107 //__________________________________________________________________________
109  TNamed(other),
110  fh3JetPtDRTrackPt(other.fh3JetPtDRTrackPt),
111  fBinWidthJetPt(other.fBinWidthJetPt),
112  fJetPtMin(other.fJetPtMin),
113  fJetPtMax(other.fJetPtMax),
114  fCollTemplates(other.fCollTemplates),
115  fInitialized(other.fInitialized),
116  fEfficiencyFixed(other.fEfficiencyFixed),
117  fhEfficiency(other.fhEfficiency),
118  fhSmoothEfficiency(other.fhSmoothEfficiency),
119  fCorrectpTtrack(other.fCorrectpTtrack),
120  fpAppliedEfficiency(other.fpAppliedEfficiency),
121  fhNmissing(other.fhNmissing),
122  fRndm(other.fRndm)
123 {
124  // Copy constructor.
125 }
126 
127 //__________________________________________________________________________
129 {
130  // Assignment
131  if (&other == this) return *this;
132  TNamed::operator=(other);
135  fJetPtMin = other.fJetPtMin;
136  fJetPtMax = other.fJetPtMax;
138  fInitialized = other.fInitialized;
140  fhEfficiency = other.fhEfficiency;
144  fhNmissing = other.fhNmissing;
145  fRndm = other.fRndm;
146  return *this;
147 }
148 
149 //__________________________________________________________________________
150 AliEmcalJet* AliEmcalJetByJetCorrection::Eval(const AliEmcalJet *jet, TClonesArray *fTracks) {
151 
152  if(!fInitialized) {
153  Printf("AliEmcalJetByJetCorrection %s not initialized",GetName());
154  return NULL;
155  }
156 
157  Int_t bin = GetJetPtBin(jet->Pt());
158  if(bin<0 || bin>fCollTemplates.GetEntriesFast()) return NULL;
159 
160  TH2D *hTemplate = static_cast<TH2D*>(fCollTemplates.At(bin));
161 
162  Double_t meanPt = GetMeanPtConstituents(jet,fTracks);
163  Double_t eff = GetEfficiency(meanPt);
164  fpAppliedEfficiency->Fill(meanPt,eff);
165 
166  Double_t fillarray[5]; //"#it{p}_{T,jet}", "N constituents added", "N_{constituents} #times (1/eff - 1)", "N_{truth}"
167  fillarray[0] = jet->Pt();
168 
169  //np is the estimation of missed tracks
170  Int_t np = TMath::FloorNint((double)jet->GetNumberOfTracks() * (1./eff -1.));
171  fillarray[2] = np;
172 
173  if(fExternalNmissed) {
174  np = fNMissedTracks; //if the number of missed tracks was calculated from external sources
175  fillarray[3] = fNMissedTracks;
176  }
177  //npc is the number of added tracks
178  Int_t npc=np; //take the particle missed as particle added
179  if(fNpPoisson){
180  npc=fRndm->Poisson(np); // smear the particle missed with a poissonian to get the number of added
181  }
182  fillarray[1] = npc;
183  fhNmissing->Fill(fillarray);
184 
185  TLorentzVector corrVec; corrVec.SetPtEtaPhiM(jet->Pt(),jet->Eta(),jet->Phi(),jet->M());
186 
187  Double_t mass = 0.13957; //pion mass
188  fArrayTrackCorr->Clear();
189  for(Int_t i = 0; i<npc; i++) {
190  Double_t r;
191  Double_t pt;
192  hTemplate->GetRandom2(r,pt);
193  Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
194  Double_t deta = r*TMath::Cos(t);
195  Double_t dphi = r*TMath::Sin(t);
196  TLorentzVector curVec;
197  curVec.SetPtEtaPhiM(pt,deta+jet->Eta(),dphi+jet->Phi(),mass);
198  new ((*fArrayTrackCorr)[i]) TLorentzVector(curVec);
199  corrVec+=curVec;
200  }
201  AliEmcalJet *jetCorr = new AliEmcalJet(corrVec.Pt(),corrVec.Eta(),corrVec.Phi(),corrVec.M());
202  Printf("ERROR: Need to implemet a fix here to set the jet acceptance correctly -> not done in the constructor AliEmcalJet(pt,eta,phi,m) just used. See JIRA https://alice.its.cern.ch/jira/browse/ALPHY-65. Use something like jet->SetJetAcceptanceType(fJetTask->FindJetAcceptanceType(eta,phi,r))");
203  return jetCorr;
204 }
205 
206 //__________________________________________________________________________
208  //Init templates
209 
210  if(!fh3JetPtDRTrackPt) {
211  Printf("%s fh3JetPtDRTrackPt not known",GetName());
212  fInitialized = kFALSE;
213  return;
214  }
215  if(!fhEfficiency && fEfficiencyFixed==1){
216  Printf("%s fhEfficiency not known",GetName());
217  fInitialized = kFALSE;
218  return;
219 
220  }
221  fRndm = new TRandom3(1234);
222  if(fhEfficiency){
223  fhSmoothEfficiency = (TH1D*) fh3JetPtDRTrackPt->ProjectionZ("fhSmoothEfficiency"); //copy the binning of the pTtrack axis
224  Int_t nBinsEf = fhEfficiency->GetXaxis()->GetNbins();
225  Int_t nBinsEfSmth = fhSmoothEfficiency->GetXaxis()->GetNbins();
226  Double_t smallBinW = 0.9;
227  Double_t pTFitRange[2]={6,100}; //reset in the loop, silent warning
228  for(Int_t ibeff=0;ibeff<nBinsEf;ibeff++){
229  Double_t bw = fhEfficiency->GetBinWidth(ibeff+1);
230  if(bw>smallBinW){
231  pTFitRange[0] = fhEfficiency->GetBinLowEdge(ibeff);
232  break;
233  }
234  }
235  pTFitRange[1] = fhEfficiency->GetBinLowEdge(nBinsEf+1);
236 
237  TF1* fitfuncEff = new TF1("fitfuncEff", "[0]+x*[1]", pTFitRange[0], pTFitRange[1]);
238  //fit function in the high pT region to smooth out fluctuations
239  fhEfficiency->Fit("fitfuncEff", "R0");
240 
241  for(Int_t i=0;i<nBinsEfSmth;i++){
242  Double_t binCentreT=fhSmoothEfficiency->GetBinCenter(i+1);
243  Int_t bin = fhEfficiency->FindBin(binCentreT);
244  Double_t binEff = fhEfficiency->GetBinContent(bin);
245 
246  //fill histogram fhSmoothEfficiency by interpolation or function
247  if(fhEfficiency->GetBinWidth(bin) > smallBinW){
248  fhSmoothEfficiency->SetBinContent(i+1,fitfuncEff->Eval(binCentreT));
249  } else {
250  Double_t effInterp = fhEfficiency->Interpolate(binCentreT);
251  fhSmoothEfficiency ->SetBinContent(i+1,effInterp);
252  fhSmoothEfficiency ->SetBinError(i+1,0.);
253 
254  }
255  }
256  }
257 
258  if(fJetPtMax>fh3JetPtDRTrackPt->GetXaxis()->GetXmax())
259  fJetPtMax = fh3JetPtDRTrackPt->GetXaxis()->GetXmax();
260 
261  if(fJetPtMin<fh3JetPtDRTrackPt->GetXaxis()->GetXmin())
262  fJetPtMin = fh3JetPtDRTrackPt->GetXaxis()->GetXmin();
263 
264  Double_t eps = 0.00001;
265  Int_t counter = 0;
267  Int_t binMin = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+eps);
268  Int_t binMax = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+fBinWidthJetPt-eps);
269 
270  fh3JetPtDRTrackPt->GetXaxis()->SetRange(binMin,binMax);
271  Int_t nBinspTtr = fh3JetPtDRTrackPt->GetZaxis()->GetNbins();
272  Int_t nBinsr = fh3JetPtDRTrackPt->GetYaxis()->GetNbins();
273  TH2D *h2 = dynamic_cast<TH2D*>(fh3JetPtDRTrackPt->Project3D("zy"));
274  if(h2){
275  h2->SetName(Form("hPtR_%.0f_%.0f",fh3JetPtDRTrackPt->GetXaxis()->GetBinLowEdge(binMin),fh3JetPtDRTrackPt->GetXaxis()->GetBinUpEdge(binMax)));
276  if(fCorrectpTtrack) {
277  //apply efficiency correction to pTtrack
278  for(Int_t ipTtr=0;ipTtr<nBinspTtr;ipTtr++){
279  for(Int_t ir=0;ir<nBinsr;ir++){
280  Double_t uncorr = h2->GetBinContent(ir+1,ipTtr+1);
281  Double_t corr = uncorr/GetEfficiency(h2->GetYaxis()->GetBinCenter(ipTtr+1));
282  h2->SetBinContent(ir+1, ipTtr+1, corr);
283  }
284 
285  }
286 
287  }
288  }
289  fCollTemplates.Add(h2);
290  fListOfOutput->Add(h2);
291  counter++;
292  }
293  // Int_t nt = TMath::FloorNint((fJetPtMax-fJetPtMin)/fBinWidthJetPt);
294  // Printf("nt: %d entries fCollTemplates: %d",nt,fCollTemplates.GetEntriesFast());
295 
296  fArrayTrackCorr = new TClonesArray("TLorentzVector", 20);
297  fInitialized = kTRUE;
298 
299 }
300 
301 //__________________________________________________________________________
303 
304  if(jetpt<fJetPtMin || jetpt>=fJetPtMax)
305  return -1;
306 
307  Int_t bin = TMath::FloorNint((jetpt - fJetPtMin)/fBinWidthJetPt);
308 
309  return bin;
310 }
311 
312 //________________________________________________________________________
314  Double_t eff = 1.;
315  if(fEfficiencyFixed<1.) return fEfficiencyFixed;
316  else if(fhSmoothEfficiency) {
317  Int_t bin = fhSmoothEfficiency->GetXaxis()->FindBin(pt);
318  eff = fhSmoothEfficiency->GetBinContent(bin);
319  //Printf("Efficiency (pT = %.2f, bin %d) = %.2f", pt, bin, eff);
320  }
321  return eff;
322 }
323 
324 //________________________________________________________________________
326 
327  if(!jet || !fTracks) return -1.;
328  if(jet->GetNumberOfTracks()<1) return -1;
329 
330  AliVParticle *vp;
331  Double_t sumPtCh = 0.;
332  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
333  vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
334  if(!vp) continue;
335  sumPtCh+=vp->Pt();
336  }
337  Double_t meanpt = sumPtCh/(double)(jet->GetNumberOfTracks());
338  return meanpt;
339 }
double Double_t
Definition: External.C:58
const char * title
Definition: MakeQAPdf.C:26
Double_t Eta() const
Definition: AliEmcalJet.h:108
TClonesArray * fArrayTrackCorr
TClonesArray containing the jet constituents after correction.
Double_t Phi() const
Definition: AliEmcalJet.h:104
AliEmcalJetByJetCorrection & operator=(const AliEmcalJetByJetCorrection &jet)
Double_t mass
AliEmcalJet * Eval(const AliEmcalJet *jet, TClonesArray *fTracks)
Bool_t fNpPoisson
draw Nmissing particle from a Poissonian with mean Nconst(1/eff-1)
TRandom * gRandom
TH1D * fhSmoothEfficiency
single particle efficiency smooth (see Init())
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
Double_t GetEfficiency(const Double_t pt) const
TProfile * fpAppliedEfficiency
Control profile efficiency.
int Int_t
Definition: External.C:63
TList * fListOfOutput
list containing all histograms
Bool_t fExternalNmissed
Set to true if want to give Nmissing from the MassStructureTask.
Definition: External.C:228
Definition: External.C:212
const Double_t ptmin
THnSparse * fhNmissing
pTjet vs number of added constituents (depends on settings) versus Nmissed constituents ...
Double_t Pt() const
Definition: AliEmcalJet.h:96
TRandom3 * fRndm
TRandom3 object.
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
Double_t fBinWidthJetPt
jet pt bin width in which to do correction
TObjArray fCollTemplates
templates (2D histos with track pT vs r)
Double_t fEfficiencyFixed
fixed efficiency for all pT and all types of tracks
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Bool_t fCorrectpTtrack
if true the templates are corrected by track efficiency
Double_t GetMeanPtConstituents(const AliEmcalJet *jet, TClonesArray *fTracks) const
Int_t GetJetPtBin(const Double_t jetpt) const
const Int_t nbins
Double_t M() const
Definition: AliEmcalJet.h:107
ClassImp(AliEmcalJetByJetCorrection) AliEmcalJetByJetCorrection
TH1 * fhEfficiency
single particle efficiency
Bool_t fInitialized
status of initialization
Int_t fNMissedTracks
Track missed in reconstruction calculated from external input (to be improved)