PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMassStructure.cxx

   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"
  35 #include "AliEmcalJetByJetCorrection.h"
  36
  37 #include "AliAODEvent.h"
  38
  39 #include "AliAnalysisTaskEmcalJetMassStructure.h"
  40
  41 ClassImp(AliAnalysisTaskEmcalJetMassStructure)
  42
  43 //________________________________________________________________________
  44 AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure() :
  45   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassStructure", kTRUE),
  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   fh3JetPtDRTrackPt(0),
  62   fpUsedEfficiency(0)
  63 {
  64   // Default constructor.
  65
  66   fh3PtDRMass                       = new TH3F*[fNcentBins];
  67   fh3PtDRRho                        = new TH3F*[fNcentBins];
  68   fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
  69   fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
  70   fh2PtMass                         = new TH2F*[fNcentBins];
  71   fh2PtMassCorr                     = new TH2F*[fNcentBins];
  72   fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];
  73
  74   for (Int_t i = 0; i < fNcentBins; i++) {
  75     fh3PtDRMass[i]       = 0;
  76     fh3PtDRRho[i]        = 0;
  77     fh3PtDRMassCorr[i]   = 0;
  78     fh3PtDRRhoCorr[i]    = 0;
  79     fh2PtMass[i]         = 0;
  80     fh2PtMassCorr[i]     = 0;
  81     fh3JetPtDRTrackPt[i] = 0;
  82   }
  83
  84   SetMakeGeneralHistograms(kTRUE);
  85 }
  86
  87 //________________________________________________________________________
  88 AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure(const char *name) :
  89   AliAnalysisTaskEmcalJet(name, kTRUE),
  90   fContainerBase(0),
  91   fMinFractionShared(0),
  92   fJetMassType(kRaw),
  93   fRandom(0),
  94   fEfficiencyFixed(1.),
  95   fCorrType(kMeanPtR),
  96   fEJetByJetCorr(0),
  97   fh3PtDRMass(0),
  98   fh3PtDRRho(0),
  99   fh3PtDRMassCorr(0),
 100   fh3PtDRRhoCorr(0),
 101   fh2PtMass(0),
 102   fh2PtMassCorr(0),
 103   fhnMassResponse(0),
 104   fhnMassResponseCorr(0),
 105   fh3JetPtDRTrackPt(0),
 106   fpUsedEfficiency(0)
 107 {
 108   // Standard constructor.
 109
 110   fh3PtDRMass                       = new TH3F*[fNcentBins];
 111   fh3PtDRRho                        = new TH3F*[fNcentBins];
 112   fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
 113   fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
 114   fh2PtMass                         = new TH2F*[fNcentBins];
 115   fh2PtMassCorr                     = new TH2F*[fNcentBins];
 116   fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];
 117
 118   for (Int_t i = 0; i < fNcentBins; i++) {
 119     fh3PtDRMass[i]       = 0;
 120     fh3PtDRRho[i]        = 0;
 121     fh3PtDRMassCorr[i]   = 0;
 122     fh3PtDRRhoCorr[i]    = 0;
 123     fh2PtMass[i]         = 0;
 124     fh2PtMassCorr[i]     = 0;
 125     fh3JetPtDRTrackPt[i] = 0;
 126   }
 127   SetMakeGeneralHistograms(kTRUE);
 128 }
 129
 130 //________________________________________________________________________
 131 AliAnalysisTaskEmcalJetMassStructure::~AliAnalysisTaskEmcalJetMassStructure()
 132 {
 133   // Destructor.
 134   delete fRandom;
 135 }
 136
 137 //________________________________________________________________________
 138 void AliAnalysisTaskEmcalJetMassStructure::UserCreateOutputObjects()
 139 {
 140   // Create user output.
 141
 142   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 143
 144   if(!fRandom) fRandom = new TRandom3(0);
 145
 146   if(fCorrType==kMeanPtR && fEJetByJetCorr) fEJetByJetCorr->Init();
 147
 148   Bool_t oldStatus = TH1::AddDirectoryStatus();
 149   TH1::AddDirectory(kFALSE);
 150
 151   const Int_t nBinsPt  = 200;
 152   const Double_t minPt = -50.;
 153   const Double_t maxPt = 150.;
 154   Double_t *binsPt = new Double_t[nBinsPt+1];
 155   for(Int_t i=0; i<=nBinsPt; i++) binsPt[i]=(Double_t)minPt + (maxPt-minPt)/nBinsPt*(Double_t)i ;
 156
 157   const Int_t nBinsM  = 120;
 158   const Double_t minM = -20.;
 159   const Double_t maxM = 40.;
 160   Double_t *binsM = new Double_t[nBinsM+1];
 161   for(Int_t i=0; i<=nBinsM; i++) binsM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;
 162
 163   const Int_t nBinsR  = 80;
 164   const Double_t minR = -0.005;
 165   const Double_t maxR = 0.795;
 166   Double_t *binsR = new Double_t[nBinsR+1];
 167   for(Int_t i=0; i<=nBinsR; i++) binsR[i]=(Double_t)minR + (maxR-minR)/nBinsR*(Double_t)i ;
 168
 169   //track pt
 170   Double_t binWidth1 = 0.1;
 171   Double_t binWidth2 = 1.;
 172   Double_t binWidth3 = 1.;
 173   const Float_t ptmin1 =  0.;
 174   const Float_t ptmax1 =  5.;
 175   const Float_t ptmin2 =  ptmax1;
 176   const Float_t ptmax2 =  10.;
 177   const Float_t ptmin3 =  ptmax2 ;
 178   const Float_t ptmax3 =  100.;
 179   const Int_t nbin11 = (int)((ptmax1-ptmin1)/binWidth1);
 180   const Int_t nbin12 = (int)((ptmax2-ptmin2)/binWidth2)+nbin11;
 181   const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
 182   Int_t nBinsPtTr=nbin13;
 183   //Create array with low edges of each bin
 184   Double_t *binsPtTr=new Double_t[nBinsPtTr+1];
 185   for(Int_t i=0; i<=nBinsPtTr; i++) {
 186     if(i<=nbin11) binsPtTr[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
 187     if(i<=nbin12 && i>nbin11) binsPtTr[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
 188     if(i<=nbin13 && i>nbin12) binsPtTr[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
 189   }
 190
 191   //Binning for THnSparse
 192   const Int_t nBinsSparse0 = 5;
 193   const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
 194   const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt};
 195   const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt};
 196
 197   TString histName = "";
 198   TString histTitle = "";
 199   for (Int_t i = 0; i < fNcentBins; i++) {
 200     histName = TString::Format("fh3PtDRMass_%d",i);
 201     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
 202     fh3PtDRMass[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 203     fOutput->Add(fh3PtDRMass[i]);
 204
 205     histName = TString::Format("fh3PtDRRho_%d",i);
 206     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
 207     fh3PtDRRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 208     fOutput->Add(fh3PtDRRho[i]);
 209
 210     histName = TString::Format("fh3PtDRMassCorr_%d",i);
 211     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
 212     fh3PtDRMassCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 213     fOutput->Add(fh3PtDRMassCorr[i]);
 214
 215     histName = TString::Format("fh3PtDRRhoCorr_%d",i);
 216     histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
 217     fh3PtDRRhoCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
 218     fOutput->Add(fh3PtDRRhoCorr[i]);
 219
 220     histName = TString::Format("fh2PtMass_%d",i);
 221     histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
 222     fh2PtMass[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 223     fOutput->Add(fh2PtMass[i]);
 224
 225     histName = TString::Format("fh2PtMassCorr_%d",i);
 226     histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
 227     fh2PtMassCorr[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
 228     fOutput->Add(fh2PtMassCorr[i]);
 229
 230     histName = TString::Format("fh3JetPtDRTrackPt_%d",i);
 231     histTitle = TString::Format("%s;#it{p}_{T,jet};r;#it{p}_{T,track}",histName.Data());
 232     fh3JetPtDRTrackPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,binsPt,nBinsR,binsR,nBinsPtTr,binsPtTr);
 233     fOutput->Add(fh3JetPtDRTrackPt[i]);
 234   }
 235
 236   histName = "fhnMassResponse";
 237   histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
 238   fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 239   fOutput->Add(fhnMassResponse);
 240
 241   histName = "fhnMassResponseCorr";
 242   histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
 243   fhnMassResponseCorr = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
 244   fOutput->Add(fhnMassResponseCorr);
 245
 246   if(fEJetByJetCorr) fpUsedEfficiency = fEJetByJetCorr->GetAppliedEfficiency();
 247   fOutput->Add(fpUsedEfficiency);
 248
 249   TH1::AddDirectory(oldStatus);
 250
 251   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 252
 253   if(binsPt)                delete [] binsPt;
 254   if(binsPtTr)              delete [] binsPtTr;
 255   if(binsM)                 delete [] binsM;
 256   if(binsR)                 delete [] binsR;
 257
 258 }
 259
 260 //________________________________________________________________________
 261 Bool_t AliAnalysisTaskEmcalJetMassStructure::Run()
 262 {
 263   // Run analysis code here, if needed. It will be executed before FillHistograms().
 264
 265   return kTRUE;
 266 }
 267
 268 //________________________________________________________________________
 269 Bool_t AliAnalysisTaskEmcalJetMassStructure::FillHistograms()
 270 {
 271   // study how jet mass builds up as function of radial distance to jet axis
 272
 273   if(fCorrType==kMeanPtR && !fEJetByJetCorr) return kFALSE;
 274
 275   AliEmcalJet* jet1 = NULL;
 276   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 277   if(jetCont) {
 278     jetCont->ResetCurrentID();
 279     while((jet1 = jetCont->GetNextAcceptJet())) {
 280
 281       Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
 282       Double_t mJet1 = GetJetMass(jet1);
 283       fh2PtMass[fCentBin]->Fill(ptJet1,mJet1);
 284       AliEmcalJet *jPart = jet1->ClosestJet();
 285       //fill detector response
 286       if(jPart) {
 287         Double_t var[4] = {mJet1,jPart->M(),ptJet1,jPart->Pt()};
 288         fhnMassResponse->Fill(var);
 289       }
 290
 291       // if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
 292       //        continue;
 293
 294       //Sort array based on distance to jet axis
 295       static Int_t indexes[9999] = {-1};
 296       static Float_t dr[9999] = {0};
 297       if(!GetSortedArray(indexes,dr,jet1)) continue;
 298
 299       for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
 300         AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
 301         if(!vp) continue;
 302         fh3JetPtDRTrackPt[fCentBin]->Fill(ptJet1,dr[indexes[i]],vp->Pt());
 303       }
 304
 305       if(fCorrType==kAnnulus) {
 306         TLorentzVector sumVec;      sumVec.SetPtEtaPhiM(0.,0.,0.,0.);
 307         TLorentzVector corrVec;     corrVec.SetPtEtaPhiM(0.,0.,0.,0.);
 308         TLorentzVector curVec;
 309         for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
 310           AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
 311           if(!vp) continue;
 312
 313           curVec.SetPtEtaPhiM(vp->Pt(),vp->Eta(),vp->Phi(),vp->M());
 314           sumVec+=curVec;
 315           corrVec+=curVec;
 316
 317           fh3PtDRMass[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.M()/mJet1);
 318           fh3PtDRRho[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.Pt()/ptJet1);
 319
 320           fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
 321           fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
 322
 323           Double_t eff = GetEfficiency(vp->Pt());
 324           Double_t rnd = fRandom->Uniform(1.);
 325           if(rnd>eff) {//put particle back at random position in annulus; using now zero width for annulus
 326             Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
 327             Double_t rr = dr[indexes[i]];
 328             rr = fRandom->Uniform(0.,jetCont->GetJetRadius());
 329             Double_t deta = rr*TMath::Cos(t);
 330             Double_t dphi = rr*TMath::Sin(t);
 331             curVec.SetPtEtaPhiM(vp->Pt(),deta+jet1->Eta(),dphi+jet1->Phi(),vp->M());
 332             corrVec+=curVec;
 333
 334             fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
 335             fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
 336           }
 337         }//track loop
 338
 339         fh2PtMassCorr[fCentBin]->Fill(corrVec.Pt(), corrVec.M());
 340         if(jPart) {
 341           Double_t varCorr[4] = {corrVec.M(),jPart->M(),corrVec.Pt(),jPart->Pt()};
 342           fhnMassResponseCorr->Fill(varCorr);
 343         }
 344       }//kAnnulus method
 345       else if(fCorrType==kMeanPtR) {
 346         //jet-by-jet correction based on templates
 347         AliEmcalJet *jetCorr = fEJetByJetCorr->Eval(jet1,jetCont->GetParticleContainer()->GetArray());
 348         if(jPart && jetCorr) {
 349           Double_t varCorr[4] = {jetCorr->M(),jPart->M(),jetCorr->Pt(),jPart->Pt()};
 350           fhnMassResponseCorr->Fill(varCorr);
 351         }
 352       }//kMeanPtR method
 353     }//jet loop
 354   }//jet container
 355
 356   return kTRUE;
 357 }
 358
 359 //________________________________________________________________________
 360 Double_t AliAnalysisTaskEmcalJetMassStructure::GetEfficiency(Double_t pt) {
 361   pt = 1.*pt;
 362   Double_t eff = 1.;
 363   if(fEfficiencyFixed<1.) return fEfficiencyFixed;
 364
 365   return eff;
 366 }
 367
 368 //________________________________________________________________________
 369 Bool_t AliAnalysisTaskEmcalJetMassStructure::GetSortedArray(Int_t indexes[], Float_t dr[], AliEmcalJet *jet) const
 370 {
 371   // sort array
 372   const Int_t n = (Int_t)jet->GetNumberOfTracks();//(Int_t)array.size();
 373   if (n < 1) return kFALSE;
 374
 375   for (Int_t i = 0; i < n; i++) {
 376     AliVParticle *vp = static_cast<AliVParticle*>(jet->TrackAt(i, fTracks));
 377     if(!vp) continue;
 378     dr[i] = jet->DeltaR(vp);
 379   }
 380   TMath::Sort(n, dr, indexes);
 381   return kTRUE;
 382 }
 383
 384
 385 //________________________________________________________________________
 386 Double_t AliAnalysisTaskEmcalJetMassStructure::GetJetMass(AliEmcalJet *jet) {
 387   //calc subtracted jet mass
 388   if(fJetMassType==kRaw)
 389     return jet->M();
 390   else if(fJetMassType==kDeriv)
 391     return jet->GetSecondOrderSubtracted();
 392
 393   return 0;
 394 }
 395
 396 //________________________________________________________________________
 397 Bool_t AliAnalysisTaskEmcalJetMassStructure::RetrieveEventObjects() {
 398   //
 399   // retrieve event objects
 400   //
 401   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 402     return kFALSE;
 403
 404   return kTRUE;
 405 }
 406
 407 //_______________________________________________________________________
 408 void AliAnalysisTaskEmcalJetMassStructure::Terminate(Option_t *)
 409 {
 410   // Called once at the end of the analysis.
 411 }
 412