2 // Utility class to do jet-by-jet correction
3 // Based on templates of ptJet vs ptTrack vs r
11 #include "TLorentzVector.h"
12 #include "AliEmcalJet.h"
13 #include "AliVParticle.h"
15 #include "AliEmcalJetByJetCorrection.h"
17 ClassImp(AliEmcalJetByJetCorrection)
19 //__________________________________________________________________________
20 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection() :
22 fh3JetPtDRTrackPt(0x0),
32 fCollTemplates.SetOwner(kTRUE);
35 //__________________________________________________________________________
36 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection(const char* name) :
38 fh3JetPtDRTrackPt(0x0),
47 // Default constructor.
48 fCollTemplates.SetOwner(kTRUE);
51 //__________________________________________________________________________
52 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection(const AliEmcalJetByJetCorrection &other) :
54 fh3JetPtDRTrackPt(other.fh3JetPtDRTrackPt),
55 fBinWidthJetPt(other.fBinWidthJetPt),
56 fJetPtMin(other.fJetPtMin),
57 fJetPtMax(other.fJetPtMax),
58 fCollTemplates(other.fCollTemplates),
59 fInitialized(other.fInitialized),
60 fEfficiencyFixed(other.fEfficiencyFixed),
61 fhEfficiency(other.fhEfficiency)
66 //__________________________________________________________________________
67 AliEmcalJetByJetCorrection& AliEmcalJetByJetCorrection::operator=(const AliEmcalJetByJetCorrection &other)
70 if (&other == this) return *this;
71 TNamed::operator=(other);
72 fh3JetPtDRTrackPt = other.fh3JetPtDRTrackPt;
73 fBinWidthJetPt = other.fBinWidthJetPt;
74 fJetPtMin = other.fJetPtMin;
75 fJetPtMax = other.fJetPtMax;
76 fCollTemplates = other.fCollTemplates;
77 fInitialized = other.fInitialized;
78 fEfficiencyFixed = other.fEfficiencyFixed;
79 fhEfficiency = other.fhEfficiency;
84 //__________________________________________________________________________
85 AliEmcalJet* AliEmcalJetByJetCorrection::Eval(const AliEmcalJet *jet, TClonesArray *fTracks) {
88 Printf("AliEmcalJetByJetCorrection %s not initialized",GetName());
92 Int_t bin = GetJetPtBin(jet->Pt());
93 if(bin<0 || bin>fCollTemplates.GetEntriesFast()) return NULL;
95 TH2D *hTemplate = static_cast<TH2D*>(fCollTemplates.At(bin));
97 Double_t meanPt = GetMeanPtConstituents(jet,fTracks);
98 Double_t eff = GetEfficiency(meanPt);
100 Int_t np = TMath::FloorNint((double)jet->GetNumberOfTracks() * (1./eff -1.));
102 TLorentzVector corrVec; corrVec.SetPtEtaPhiM(jet->Pt(),jet->Eta(),jet->Phi(),jet->M());
104 Double_t mass = 0.13957; //pion mass
106 for(Int_t i = 0; i<np; i++) {
109 hTemplate->GetRandom2(r,pt);
110 Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
111 Double_t deta = r*TMath::Cos(t);
112 Double_t dphi = r*TMath::Sin(t);
113 TLorentzVector curVec;
114 curVec.SetPtEtaPhiM(pt,deta+jet->Eta(),dphi+jet->Phi(),mass);
118 AliEmcalJet *jetCorr = new AliEmcalJet(corrVec.Pt(),corrVec.Eta(),corrVec.Phi(),corrVec.M());
123 //__________________________________________________________________________
124 void AliEmcalJetByJetCorrection::Init() {
127 if(!fh3JetPtDRTrackPt) {
128 Printf("%s fh3JetPtDRTrackPt not known",GetName());
129 fInitialized = kFALSE;
133 if(fJetPtMax>fh3JetPtDRTrackPt->GetXaxis()->GetXmax())
134 fJetPtMax = fh3JetPtDRTrackPt->GetXaxis()->GetXmax();
136 if(fJetPtMin<fh3JetPtDRTrackPt->GetXaxis()->GetXmin())
137 fJetPtMin = fh3JetPtDRTrackPt->GetXaxis()->GetXmin();
139 Double_t eps = 0.00001;
141 for(Double_t ptmin = fJetPtMin; ptmin<fJetPtMax; ptmin+=fBinWidthJetPt) {
142 Int_t binMin = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+eps);
143 Int_t binMax = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+fBinWidthJetPt-eps);
144 // Printf("%d bins: %d - %d -> %f - %f",counter,binMin,binMax,fh3JetPtDRTrackPt->GetXaxis()->GetBinLowEdge(binMin),fh3JetPtDRTrackPt->GetXaxis()->GetBinUpEdge(binMax));
146 fh3JetPtDRTrackPt->GetXaxis()->SetRange(binMin,binMax);
147 TH2D *h2 = dynamic_cast<TH2D*>(fh3JetPtDRTrackPt->Project3D("zy"));
148 h2->SetName(Form("hPtR_%.0f_%.0f",fh3JetPtDRTrackPt->GetXaxis()->GetBinLowEdge(binMin),fh3JetPtDRTrackPt->GetXaxis()->GetBinUpEdge(binMax)));
149 fCollTemplates.Add(h2);
152 // Int_t nt = TMath::FloorNint((fJetPtMax-fJetPtMin)/fBinWidthJetPt);
153 // Printf("nt: %d entries fCollTemplates: %d",nt,fCollTemplates.GetEntriesFast());
155 fInitialized = kTRUE;
159 //__________________________________________________________________________
160 Int_t AliEmcalJetByJetCorrection::GetJetPtBin(const Double_t jetpt) const {
162 if(jetpt<fJetPtMin || jetpt>=fJetPtMax)
165 Int_t bin = TMath::FloorNint((jetpt - fJetPtMin)/fBinWidthJetPt);
170 //________________________________________________________________________
171 Double_t AliEmcalJetByJetCorrection::GetEfficiency(const Double_t pt) const {
173 if(fEfficiencyFixed<1.) return fEfficiencyFixed;
174 else if(fhEfficiency) {
175 Int_t bin = fhEfficiency->GetXaxis()->FindBin(pt);
176 eff = fhEfficiency->GetBinContent(bin);
181 //________________________________________________________________________
182 Double_t AliEmcalJetByJetCorrection::GetMeanPtConstituents(const AliEmcalJet *jet, TClonesArray *fTracks) const {
184 if(!jet || !fTracks) return -1.;
185 if(jet->GetNumberOfTracks()<1) return -1;
188 Double_t sumPtCh = 0.;
189 for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
190 vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
194 Double_t meanpt = sumPtCh/(double)(jet->GetNumberOfTracks());