[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetMassStructure.cxx

//
// Jet mass structure analysis task.
//
// Author: M.Verweij

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TProfile2D.h>
#include <THnSparse.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TProfile.h>
#include <TChain.h>
#include <TSystem.h>
#include <TFile.h>
#include <TKey.h>
#include <TMath.h>
#include <TRandom3.h>

#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliEmcalParticle.h"
#include "AliMCEvent.h"
#include "AliGenPythiaEventHeader.h"
#include "AliAODMCHeader.h"
#include "AliMCEvent.h"
#include "AliAnalysisManager.h"
#include "AliJetContainer.h"

#include "AliAODEvent.h"

#include "AliAnalysisTaskEmcalJetMassStructure.h"

ClassImp(AliAnalysisTaskEmcalJetMassStructure)

//________________________________________________________________________
AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassStructure", kTRUE),
  fContainerBase(0),
  fMinFractionShared(0),
  fJetMassType(kRaw),
  fRandom(0),
  fEfficiencyFixed(1.),
  fh3PtDRMass(0),
  fh3PtDRRho(0),
  fh3PtDRMassCorr(0),
  fh3PtDRRhoCorr(0),
  fh2PtMass(0),
  fh2PtMassCorr(0),
  fhnMassResponse(0),
  fhnMassResponseCorr(0),
  fh3JetPtDRTrackPt(0)
{
  // Default constructor.

  fh3PtDRMass                       = new TH3F*[fNcentBins];
  fh3PtDRRho                        = new TH3F*[fNcentBins];
  fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
  fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
  fh2PtMass                         = new TH2F*[fNcentBins];
  fh2PtMassCorr                     = new TH2F*[fNcentBins];
  fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtDRMass[i]       = 0;
    fh3PtDRRho[i]        = 0;
    fh3PtDRMassCorr[i]   = 0;
    fh3PtDRRhoCorr[i]    = 0;
    fh2PtMass[i]         = 0;
    fh2PtMassCorr[i]     = 0;
    fh3JetPtDRTrackPt[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetMassStructure::AliAnalysisTaskEmcalJetMassStructure(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),  
  fContainerBase(0),
  fMinFractionShared(0),
  fJetMassType(kRaw),
  fRandom(0),
  fEfficiencyFixed(1.),
  fh3PtDRMass(0),
  fh3PtDRRho(0),
  fh3PtDRMassCorr(0),
  fh3PtDRRhoCorr(0),
  fh2PtMass(0),
  fh2PtMassCorr(0),
  fhnMassResponse(0),
  fhnMassResponseCorr(0),
  fh3JetPtDRTrackPt(0)
{
  // Standard constructor.

  fh3PtDRMass                       = new TH3F*[fNcentBins];
  fh3PtDRRho                        = new TH3F*[fNcentBins];
  fh3PtDRMassCorr                   = new TH3F*[fNcentBins];
  fh3PtDRRhoCorr                    = new TH3F*[fNcentBins];
  fh2PtMass                         = new TH2F*[fNcentBins];
  fh2PtMassCorr                     = new TH2F*[fNcentBins];
  fh3JetPtDRTrackPt                 = new TH3F*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fh3PtDRMass[i]       = 0;
    fh3PtDRRho[i]        = 0; 
    fh3PtDRMassCorr[i]   = 0;
    fh3PtDRRhoCorr[i]    = 0;
    fh2PtMass[i]         = 0;
    fh2PtMassCorr[i]     = 0;
    fh3JetPtDRTrackPt[i] = 0;
  }

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetMassStructure::~AliAnalysisTaskEmcalJetMassStructure()
{
  // Destructor.

  delete fRandom;
}

//________________________________________________________________________
void AliAnalysisTaskEmcalJetMassStructure::UserCreateOutputObjects()
{
  // Create user output.

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  if(!fRandom) fRandom = new TRandom3(0);

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);

  const Int_t nBinsPt  = 200;
  const Double_t minPt = -50.;
  const Double_t maxPt = 150.;
  Double_t *binsPt = new Double_t[nBinsPt+1];
  for(Int_t i=0; i<=nBinsPt; i++) binsPt[i]=(Double_t)minPt + (maxPt-minPt)/nBinsPt*(Double_t)i ;

  const Int_t nBinsM  = 120;
  const Double_t minM = -20.;
  const Double_t maxM = 40.;
  Double_t *binsM = new Double_t[nBinsM+1];
  for(Int_t i=0; i<=nBinsM; i++) binsM[i]=(Double_t)minM + (maxM-minM)/nBinsM*(Double_t)i ;

  const Int_t nBinsR  = 80;
  const Double_t minR = -0.005;
  const Double_t maxR = 0.795;
  Double_t *binsR = new Double_t[nBinsR+1];
  for(Int_t i=0; i<=nBinsR; i++) binsR[i]=(Double_t)minR + (maxR-minR)/nBinsR*(Double_t)i ;

  //track pt
  Double_t binWidth1 = 0.1;
  Double_t binWidth2 = 1.;
  Double_t binWidth3 = 1.;
  const Float_t ptmin1 =  0.;
  const Float_t ptmax1 =  5.;
  const Float_t ptmin2 =  ptmax1;
  const Float_t ptmax2 =  10.;
  const Float_t ptmin3 =  ptmax2 ;
  const Float_t ptmax3 =  100.;
  const Int_t nbin11 = (int)((ptmax1-ptmin1)/binWidth1);
  const Int_t nbin12 = (int)((ptmax2-ptmin2)/binWidth2)+nbin11;
  const Int_t nbin13 = (int)((ptmax3-ptmin3)/binWidth3)+nbin12;
  Int_t nBinsPtTr=nbin13;
  //Create array with low edges of each bin
  Double_t *binsPtTr=new Double_t[nBinsPtTr+1];
  for(Int_t i=0; i<=nBinsPtTr; i++) {
    if(i<=nbin11) binsPtTr[i]=(Double_t)ptmin1 + (ptmax1-ptmin1)/nbin11*(Double_t)i ;
    if(i<=nbin12 && i>nbin11) binsPtTr[i]=(Double_t)ptmin2 + (ptmax2-ptmin2)/(nbin12-nbin11)*((Double_t)i-(Double_t)nbin11) ;
    if(i<=nbin13 && i>nbin12) binsPtTr[i]=(Double_t)ptmin3 + (ptmax3-ptmin3)/(nbin13-nbin12)*((Double_t)i-(Double_t)nbin12) ;
  }

  //Binning for THnSparse
  const Int_t nBinsSparse0 = 5;
  const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
  const Double_t xmin0[nBinsSparse0]  = { minM, minM, minPt, minPt};
  const Double_t xmax0[nBinsSparse0]  = { maxM, maxM, maxPt, maxPt};

  TString histName = "";
  TString histTitle = "";
  for (Int_t i = 0; i < fNcentBins; i++) {
    histName = TString::Format("fh3PtDRMass_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
    fh3PtDRMass[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
    fOutput->Add(fh3PtDRMass[i]);

    histName = TString::Format("fh3PtDRRho_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
    fh3PtDRRho[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
    fOutput->Add(fh3PtDRRho[i]);

    histName = TString::Format("fh3PtDRMassCorr_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum m",histName.Data());
    fh3PtDRMassCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
    fOutput->Add(fh3PtDRMassCorr[i]);

    histName = TString::Format("fh3PtDRRhoCorr_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};r;sum #it{p}_{T}",histName.Data());
    fh3PtDRRhoCorr[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsR,minR,maxR,101,-0.005,1.005);
    fOutput->Add(fh3PtDRRhoCorr[i]);

    histName = TString::Format("fh2PtMass_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
    fh2PtMass[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
    fOutput->Add(fh2PtMass[i]);

    histName = TString::Format("fh2PtMassCorr_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};#it{M}_{jet}",histName.Data());
    fh2PtMassCorr[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsM,minM,maxM);
    fOutput->Add(fh2PtMassCorr[i]);

    histName = TString::Format("fh3JetPtDRTrackPt_%d",i);
    histTitle = TString::Format("%s;#it{p}_{T,jet};r;#it{p}_{T,track}",histName.Data());
    fh3JetPtDRTrackPt[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,binsPt,nBinsR,binsR,nBinsPtTr,binsPtTr);
    fOutput->Add(fh3JetPtDRTrackPt[i]);
  }

  histName = "fhnMassResponse";
  histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
  fhnMassResponse = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
  fOutput->Add(fhnMassResponse);

  histName = "fhnMassResponseCorr";
  histTitle = Form("%s;#it{M}_{det};#it{M}_{part};#it{p}_{T,det};#it{p}_{T,part}",histName.Data());
  fhnMassResponseCorr = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
  fOutput->Add(fhnMassResponseCorr);

  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
    if (h1){
      h1->Sumw2();
      continue;
    }
    THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
    if(hn)hn->Sumw2();
  }

  TH1::AddDirectory(oldStatus);

  PostData(1, fOutput); // Post data for ALL output slots > 0 here.

  if(binsPt)                delete [] binsPt;
  if(binsPtTr)              delete [] binsPtTr;
  if(binsM)                 delete [] binsM;
  if(binsR)                 delete [] binsR;

}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetMassStructure::Run()
{
  // Run analysis code here, if needed. It will be executed before FillHistograms().

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetMassStructure::FillHistograms()
{
  // Fill histograms.

  // study how jet mass builds up as function of radial distance to jet axis

  AliEmcalJet* jet1 = NULL;
  AliJetContainer *jetCont = GetJetContainer(fContainerBase);
  if(jetCont) {
    jetCont->ResetCurrentID();
    while((jet1 = jetCont->GetNextAcceptJet())) {

      Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
      Double_t mJet1 = GetJetMass(jet1);
           
      // if(jet1->GetTagStatus()<1 || !jet1->GetTaggedJet())
      // 	continue;
     
      //Sort array based on distance to jet axis
      static Int_t indexes[9999] = {-1};
      static Float_t dr[9999] = {0};
      if(!GetSortedArray(indexes,dr,jet1)) continue;
 
      AliVParticle *vp;
      TLorentzVector sumVec;      sumVec.SetPtEtaPhiM(0.,0.,0.,0.);
      TLorentzVector corrVec;     corrVec.SetPtEtaPhiM(0.,0.,0.,0.);
      TLorentzVector curVec;
      for(Int_t i=jet1->GetNumberOfTracks()-1; i>-1; i--) {
	vp = static_cast<AliVParticle*>(jet1->TrackAt(indexes[i], fTracks));
	if(!vp) continue;

	fh3JetPtDRTrackPt[fCentBin]->Fill(ptJet1,dr[indexes[i]],vp->Pt());

	curVec.SetPtEtaPhiM(vp->Pt(),vp->Eta(),vp->Phi(),vp->M());
	sumVec+=curVec;
	corrVec+=curVec;
	//	Printf("%d  %f",i,dr[indexes[i]]);
	fh3PtDRMass[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.M()/mJet1);
	fh3PtDRRho[fCentBin]->Fill(ptJet1,dr[indexes[i]],sumVec.Pt()/ptJet1);

	fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
	fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);

	Double_t eff = GetEfficiency(vp->Pt());
	Double_t rnd = fRandom->Uniform(1.);
	if(rnd>eff) {//put particle back at random position in annulus; using now zero width for annulus
	  Double_t t = TMath::TwoPi()*gRandom->Uniform(1.);
	  Double_t rr = dr[indexes[i]];
	  rr = fRandom->Uniform(0.,jetCont->GetJetRadius());
	  Double_t deta = rr*TMath::Cos(t);
	  Double_t dphi = rr*TMath::Sin(t);
	  curVec.SetPtEtaPhiM(vp->Pt(),deta+jet1->Eta(),dphi+jet1->Phi(),vp->M());
	  corrVec+=curVec;

	  fh3PtDRMassCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.M()/mJet1);
	  fh3PtDRRhoCorr[fCentBin]->Fill(ptJet1,dr[indexes[i]],corrVec.Pt()/ptJet1);
	}
      }

      fh2PtMass[fCentBin]->Fill(ptJet1,mJet1);
      fh2PtMassCorr[fCentBin]->Fill(corrVec.Pt(), corrVec.M());

      //fill detector response
      AliEmcalJet *jPart = jet1->ClosestJet();
      if(jPart) {
	Double_t var[4] = {mJet1,jPart->M(),ptJet1,jPart->Pt()};
	fhnMassResponse->Fill(var);

	Double_t varCorr[4] = {corrVec.M(),jPart->M(),corrVec.Pt(),jPart->Pt()};
	fhnMassResponseCorr->Fill(varCorr);
      }
    }
  }
  
  return kTRUE;
}

//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetMassStructure::GetEfficiency(Double_t pt) {
  pt = 1.*pt;
  Double_t eff = 1.;
  if(fEfficiencyFixed<1.) return fEfficiencyFixed;

  return eff;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetMassStructure::GetSortedArray(Int_t indexes[], Float_t dr[], AliEmcalJet *jet) const
{
  // sort array

  const Int_t n = (Int_t)jet->GetNumberOfTracks();//(Int_t)array.size();

  if (n < 1)
    return kFALSE;

  AliVParticle *vp;
  for (Int_t i = 0; i < n; i++) {
    vp = static_cast<AliVParticle*>(jet->TrackAt(i, fTracks));
    if(!vp) continue;
    dr[i] = jet->DeltaR(vp);
  }
   
  TMath::Sort(n, dr, indexes);

  return kTRUE;
}


//________________________________________________________________________
Double_t AliAnalysisTaskEmcalJetMassStructure::GetJetMass(AliEmcalJet *jet) {
  //calc subtracted jet mass
  if(fJetMassType==kRaw)
    return jet->M();
  else if(fJetMassType==kDeriv)
    return jet->GetSecondOrderSubtracted();
  
  return 0;
}

//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetMassStructure::RetrieveEventObjects() {
  //
  // retrieve event objects
  //
  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
    return kFALSE;

  return kTRUE;
}

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