PWGJE/EMCALJetTasks/AliEmcalJetByJetCorrection.cxx

   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
  16 #include "AliEmcalJetByJetCorrection.h"
  17
  18 ClassImp(AliEmcalJetByJetCorrection)
  19
  20 //__________________________________________________________________________
  21 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection() :
  22 TNamed(),
  23   fh3JetPtDRTrackPt(0x0),
  24   fBinWidthJetPt(10.),
  25   fJetPtMin(0.),
  26   fJetPtMax(150.),
  27   fCollTemplates(),
  28   fInitialized(kFALSE),
  29   fEfficiencyFixed(1.),
  30   fhEfficiency(0),
  31   fpAppliedEfficiency(0)
  32 {
  33   // Dummy constructor.
  34   fCollTemplates.SetOwner(kTRUE);
  35 }
  36
  37 //__________________________________________________________________________
  38 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection(const char* name) :
  39   TNamed(name, name),
  40   fh3JetPtDRTrackPt(0x0),
  41   fBinWidthJetPt(10.),
  42   fJetPtMin(0.),
  43   fJetPtMax(150.),
  44   fCollTemplates(),
  45   fInitialized(kFALSE),
  46   fEfficiencyFixed(1.),
  47   fhEfficiency(0),
  48   fpAppliedEfficiency(0)
  49 {
  50   // Default constructor.
  51   fCollTemplates.SetOwner(kTRUE);
  52
  53   const Int_t nBinPt = 118; //0-2: 20 bins; 2-100: 98 bins
  54   Double_t binLimitsPt[nBinPt+1];
  55   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
  56     if(iPt<20){
  57       binLimitsPt[iPt] = 0. + (Double_t)iPt*0.15;
  58     } else {// 1.0
  59       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 1.0;
  60     }
  61   }
  62   fpAppliedEfficiency = new TProfile("fpAppliedEfficiency","fpAppliedEfficiency",nBinPt,binLimitsPt);
  63 }
  64
  65 //__________________________________________________________________________
  66 AliEmcalJetByJetCorrection::AliEmcalJetByJetCorrection(const AliEmcalJetByJetCorrection &other) :
  67   TNamed(other),
  68   fh3JetPtDRTrackPt(other.fh3JetPtDRTrackPt),
  69   fBinWidthJetPt(other.fBinWidthJetPt),
  70   fJetPtMin(other.fJetPtMin),
  71   fJetPtMax(other.fJetPtMax),
  72   fCollTemplates(other.fCollTemplates),
  73   fInitialized(other.fInitialized),
  74   fEfficiencyFixed(other.fEfficiencyFixed),
  75   fhEfficiency(other.fhEfficiency),
  76   fpAppliedEfficiency(other.fpAppliedEfficiency)
  77 {
  78   // Copy constructor.
  79 }
  80
  81 //__________________________________________________________________________
  82 AliEmcalJetByJetCorrection& AliEmcalJetByJetCorrection::operator=(const AliEmcalJetByJetCorrection &other)
  83 {
  84   // Assignment
  85   if (&other == this) return *this;
  86   TNamed::operator=(other);
  87   fh3JetPtDRTrackPt = other.fh3JetPtDRTrackPt;
  88   fBinWidthJetPt     = other.fBinWidthJetPt;
  89   fJetPtMin          = other.fJetPtMin;
  90   fJetPtMax          = other.fJetPtMax;
  91   fCollTemplates     = other.fCollTemplates;
  92   fInitialized       = other.fInitialized;
  93   fEfficiencyFixed   = other.fEfficiencyFixed;
  94   fhEfficiency       = other.fhEfficiency;
  95   fpAppliedEfficiency= other.fpAppliedEfficiency;
  96
  97   return *this;
  98 }
  99
 100 //__________________________________________________________________________
 101 AliEmcalJet* AliEmcalJetByJetCorrection::Eval(const AliEmcalJet *jet, TClonesArray *fTracks) {
 102
 103   if(!fInitialized) {
 104     Printf("AliEmcalJetByJetCorrection %s not initialized",GetName());
 105     return NULL;
 106   }
 107
 108   Int_t bin = GetJetPtBin(jet->Pt());
 109   if(bin<0 || bin>fCollTemplates.GetEntriesFast()) return NULL;
 110
 111   TH2D *hTemplate = static_cast<TH2D*>(fCollTemplates.At(bin));
 112
 113   Double_t meanPt = GetMeanPtConstituents(jet,fTracks);
 114   Double_t eff = GetEfficiency(meanPt);
 115   fpAppliedEfficiency->Fill(meanPt,eff);
 116
 117   Int_t np = TMath::FloorNint((double)jet->GetNumberOfTracks() * (1./eff -1.));
 118
 119   TLorentzVector corrVec; corrVec.SetPtEtaPhiM(jet->Pt(),jet->Eta(),jet->Phi(),jet->M());
 120
 121   Double_t mass = 0.13957; //pion mass
 122
 123   for(Int_t i = 0; i<np; i++) {
 124     Double_t r;
 125     Double_t pt;
 126     hTemplate->GetRandom2(r,pt);
 127     Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
 128     Double_t deta = r*TMath::Cos(t);
 129     Double_t dphi = r*TMath::Sin(t);
 130     TLorentzVector curVec;
 131     curVec.SetPtEtaPhiM(pt,deta+jet->Eta(),dphi+jet->Phi(),mass);
 132     corrVec+=curVec;
 133   }
 134
 135   AliEmcalJet *jetCorr = new AliEmcalJet(corrVec.Pt(),corrVec.Eta(),corrVec.Phi(),corrVec.M());
 136
 137   return jetCorr;
 138 }
 139
 140 //__________________________________________________________________________
 141 void AliEmcalJetByJetCorrection::Init() {
 142   //Init templates
 143
 144   if(!fh3JetPtDRTrackPt) {
 145     Printf("%s fh3JetPtDRTrackPt not known",GetName());
 146     fInitialized = kFALSE;
 147     return;
 148   }
 149
 150   if(fJetPtMax>fh3JetPtDRTrackPt->GetXaxis()->GetXmax())
 151     fJetPtMax = fh3JetPtDRTrackPt->GetXaxis()->GetXmax();
 152
 153   if(fJetPtMin<fh3JetPtDRTrackPt->GetXaxis()->GetXmin())
 154     fJetPtMin = fh3JetPtDRTrackPt->GetXaxis()->GetXmin();
 155
 156   Double_t eps = 0.00001;
 157   Int_t counter = 0;
 158   for(Double_t ptmin = fJetPtMin; ptmin<fJetPtMax; ptmin+=fBinWidthJetPt) {
 159     Int_t binMin = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+eps);
 160     Int_t binMax = fh3JetPtDRTrackPt->GetXaxis()->FindBin(ptmin+fBinWidthJetPt-eps);
 161     //    Printf("%d bins: %d - %d -> %f - %f",counter,binMin,binMax,fh3JetPtDRTrackPt->GetXaxis()->GetBinLowEdge(binMin),fh3JetPtDRTrackPt->GetXaxis()->GetBinUpEdge(binMax));
 162
 163     fh3JetPtDRTrackPt->GetXaxis()->SetRange(binMin,binMax);
 164     TH2D *h2 = dynamic_cast<TH2D*>(fh3JetPtDRTrackPt->Project3D("zy"));
 165     h2->SetName(Form("hPtR_%.0f_%.0f",fh3JetPtDRTrackPt->GetXaxis()->GetBinLowEdge(binMin),fh3JetPtDRTrackPt->GetXaxis()->GetBinUpEdge(binMax)));
 166     fCollTemplates.Add(h2);
 167     counter++;
 168   }
 169   //  Int_t nt = TMath::FloorNint((fJetPtMax-fJetPtMin)/fBinWidthJetPt);
 170   //  Printf("nt: %d entries fCollTemplates: %d",nt,fCollTemplates.GetEntriesFast());
 171
 172   fInitialized = kTRUE;
 173
 174 }
 175
 176 //__________________________________________________________________________
 177 Int_t AliEmcalJetByJetCorrection::GetJetPtBin(const Double_t jetpt) const {
 178
 179   if(jetpt<fJetPtMin || jetpt>=fJetPtMax)
 180     return -1;
 181
 182   Int_t bin = TMath::FloorNint((jetpt - fJetPtMin)/fBinWidthJetPt);
 183
 184   return bin;
 185 }
 186
 187 //________________________________________________________________________
 188 Double_t AliEmcalJetByJetCorrection::GetEfficiency(const Double_t pt) const {
 189   Double_t eff = 1.;
 190   if(fEfficiencyFixed<1.) return fEfficiencyFixed;
 191   else if(fhEfficiency) {
 192     Int_t bin = fhEfficiency->GetXaxis()->FindBin(pt);
 193     eff = fhEfficiency->GetBinContent(bin);
 194   }
 195   return eff;
 196 }
 197
 198 //________________________________________________________________________
 199 Double_t AliEmcalJetByJetCorrection::GetMeanPtConstituents(const AliEmcalJet *jet, TClonesArray *fTracks) const {
 200
 201   if(!jet || !fTracks) return -1.;
 202   if(jet->GetNumberOfTracks()<1) return -1;
 203
 204   AliVParticle *vp;
 205   Double_t sumPtCh = 0.;
 206   for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 207     vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 208     if(!vp) continue;
 209     sumPtCh+=vp->Pt();
 210   }
 211   Double_t meanpt = sumPtCh/(double)(jet->GetNumberOfTracks());
 212   return meanpt;
 213 }