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

// $Id$
//
// Dcal dijet performance task
//
// Author: R. Reed

#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THnSparse.h>
#include <TList.h>
#include <TLorentzVector.h>

#include "AliVCluster.h"
#include "AliAODCaloCluster.h"
#include "AliESDCaloCluster.h"
#include "AliVTrack.h"
#include "AliEmcalJet.h"
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliJetContainer.h"
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"
#include "AliPicoTrack.h"

#include "AliAnalysisTaskDcalDijetPerf.h"

ClassImp(AliAnalysisTaskDcalDijetPerf)

//________________________________________________________________________
AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf() : 
  AliAnalysisTaskEmcalJet("AliAnalysisTaskDcalDijetPerf", kTRUE),
  fHistTracksPt(0),
  fHistTracksEtaPhi(0),
  fHistClustersPt(0),
  fHistLeadingJetPt(0),
  fHistJetsPhiEta(0),
  fHistJetsPtArea(0),
  fHistJetsPtLeadHad(0),
  fHistJetsCorrPtArea(0),
  fHistJet1(0),
  fHistJet1m(0),
  fHistJet1nm(0),
  fHistJet2(0),
  fHistJet1to2(0),
  fJetsCont(0),
  fJetsCont2(0),
  fTracksCont(0),
  fCaloClustersCont(0)

{
  // Default constructor.

  fHistTracksPt       = new TH1*[fNcentBins];
  fHistTracksEtaPhi   = new TH2*[fNcentBins];
  fHistClustersPt     = new TH1*[fNcentBins];
  fHistLeadingJetPt   = new TH1*[fNcentBins];
  fHistJetsPhiEta     = new TH2*[fNcentBins];
  fHistJetsPtArea     = new TH2*[fNcentBins];
  fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  fHistJetsCorrPtArea = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistTracksPt[i] = 0;
    fHistTracksEtaPhi[i] = 0;
    fHistClustersPt[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHistJetsPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtLeadHad[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
  }
    fHistJet1 = 0;
    fHistJet1m = 0;
    fHistJet1nm = 0;
    fHistJet2 = 0;
    fHistJet1to2 = 0;
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf(const char *name) : 
  AliAnalysisTaskEmcalJet(name, kTRUE),
  fHistTracksPt(0),
  fHistTracksEtaPhi(0),
  fHistClustersPt(0),
  fHistLeadingJetPt(0),
  fHistJetsPhiEta(0),
  fHistJetsPtArea(0),
  fHistJetsPtLeadHad(0),
  fHistJetsCorrPtArea(0),
  fHistJet1(0),
  fHistJet1m(0),
  fHistJet1nm(0),
  fHistJet2(0),
  fHistJet1to2(0),
  fJetsCont(0),
  fJetsCont2(0),
  fTracksCont(0),
  fCaloClustersCont(0)
{
  // Standard constructor.

  fHistTracksPt       = new TH1*[fNcentBins];
  fHistTracksEtaPhi   = new TH2*[fNcentBins];
  fHistClustersPt     = new TH1*[fNcentBins];
  fHistLeadingJetPt   = new TH1*[fNcentBins];
  fHistJetsPhiEta     = new TH2*[fNcentBins];
  fHistJetsPtArea     = new TH2*[fNcentBins];
  fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  fHistJetsCorrPtArea = new TH2*[fNcentBins];

  for (Int_t i = 0; i < fNcentBins; i++) {
    fHistTracksPt[i] = 0;
    fHistTracksEtaPhi[i] = 0;
    fHistClustersPt[i] = 0;
    fHistLeadingJetPt[i] = 0;
    fHistJetsPhiEta[i] = 0;
    fHistJetsPtArea[i] = 0;
    fHistJetsPtLeadHad[i] = 0;
    fHistJetsCorrPtArea[i] = 0;
  }
  fHistJet1 = 0;
  fHistJet1m = 0;
  fHistJet1nm = 0;
  fHistJet2 = 0;
  fHistJet1to2 = 0;

  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskDcalDijetPerf::~AliAnalysisTaskDcalDijetPerf()
{
  // Destructor.
}

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

  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();

  fJetsCont           = GetJetContainer(0);
  fJetsCont2           = GetJetContainer(1);
  if(fJetsCont) { //get particles and clusters connected to jets
    fTracksCont       = fJetsCont->GetParticleContainer();
    fCaloClustersCont = fJetsCont->GetClusterContainer();
  } else {        //no jets, just analysis tracks and clusters
    fTracksCont       = GetParticleContainer(0);
    fCaloClustersCont = GetClusterContainer(0);
  }
  fTracksCont->SetClassName("AliVTrack");
  //fCaloClustersCont->SetClassName("AliVCluster");
  TString histname;

  for (Int_t i = 0; i < fNcentBins; i++) {
    if (fParticleCollArray.GetEntriesFast()>0) {
      histname = "fHistTracksPt_";
      histname += i;
      fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
      fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistTracksPt[i]);
      histname = "fHistTracksEtaPhi_";
      histname += i;
        fHistTracksEtaPhi[i] = new TH2F(histname.Data(), histname.Data(), fNbins, -0.7, 0.7, fNbins, 0, TMath::Pi()*2);
      fHistTracksEtaPhi[i]->GetXaxis()->SetTitle("eta");
      fHistTracksEtaPhi[i]->GetYaxis()->SetTitle("phi");
      fOutput->Add(fHistTracksEtaPhi[i]);
    }

    if (fClusterCollArray.GetEntriesFast()>0) {
      histname = "fHistClustersPt_";
      histname += i;
      fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
      fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistClustersPt[i]);
    }

    if (fJetCollArray.GetEntriesFast()>0) {
      histname = "fHistLeadingJetPt_";
      histname += i;
      fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
      fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
      fOutput->Add(fHistLeadingJetPt[i]);
      
      histname = "fHistJetsPhiEta_";
      histname += i;
      fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
      fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
      fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
      fOutput->Add(fHistJetsPhiEta[i]);
      
      histname = "fHistJetsPtArea_";
      histname += i;
      fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
      fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
      fOutput->Add(fHistJetsPtArea[i]);

      histname = "fHistJetsPtLeadHad_";
      histname += i;
      fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
      fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
      fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
      fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
      fOutput->Add(fHistJetsPtLeadHad[i]);
    
      if (!(GetJetContainer()->GetRhoName().IsNull())) {
	histname = "fHistJetsCorrPtArea_";
	histname += i;
	fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
	fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
	fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
	fOutput->Add(fHistJetsCorrPtArea[i]);
      }
    }
  }

  Int_t nbins[] = {150,100,100,100,150};
  Double_t xmin[] = {  0,-0.7,             0, 0,   0};
  Double_t xmax[] = {150, 0.7,TMath::TwoPi(), 1, 150};
  fHistJet1 = new THnSparseF("Jets1Collection","Jets1Collection",5,nbins,xmin,xmax);
  fOutput->Add(fHistJet1);
  fHistJet1m = new THnSparseF("Jets1CollectionMatched","Jets1Collection",5,nbins,xmin,xmax);
  fOutput->Add(fHistJet1m);
  fHistJet1nm = new THnSparseF("Jets1CollectionNotMatched","Jets1Collection",5,nbins,xmin,xmax);
  fOutput->Add(fHistJet1nm);
  fHistJet2 = new THnSparseF("Jets2Collection","Jets2Collection",5,nbins,xmin,xmax);
  fOutput->Add(fHistJet2);
  
  Int_t nbins2[] = {150,100,100,100,150,100,100,100,100};
  Double_t xmin2[] = {0,-0.7,0,0,0,-0.7,0,0,0};
  Double_t xmax2[] = {150,0.7,6.28,1,150,0.7,6.28,1,0.2};
  fHistJet1to2 = new THnSparseF("Jets1to2Collection","Jets1to2Collection",9,nbins2,xmin2,xmax2);
  fOutput->Add(fHistJet1to2);
  
  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
}

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

  if (fTracksCont) {
    AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
    while(track) {
      fHistTracksPt[fCentBin]->Fill(track->Pt());
      fHistTracksEtaPhi[fCentBin]->Fill(track->Eta(),track->Phi());
      track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
    }
  }
  
  if (fCaloClustersCont) {
    AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
    while(cluster) {
      TLorentzVector nPart;
      cluster->GetMomentum(nPart, fVertex);
      fHistClustersPt[fCentBin]->Fill(nPart.Pt());

      cluster = fCaloClustersCont->GetNextAcceptCluster();
    }
  }

    int N1 = 0;
  if (fJetsCont) {
    AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0); 
    while(jet) {
        Float_t NEFpT = jet->Pt()*jet->NEF();
        N1++;
        //cout<<"loop 1 jet 1 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<" and NEFpT "<<NEFpT<<endl;
      fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
      fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
      Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
      fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
        
        //110 degrees in the azimuthal angle |eta|<0.7
        //EMCal Tower 0.0143 x 0.0143
        //16x16 =  0.2288 x 0.2288 -> R = 0.1144
        
        Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
        //cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
        fHistJet1->Fill(jetarray);
      if (fHistJetsCorrPtArea[fCentBin]) {
	Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
	fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
      }
      jet = fJetsCont->GetNextAcceptJet();
    }
    
    jet = fJetsCont->GetLeadingJet();
    if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
  }
   // cout<<"DONE LOOP 1"<<endl;
    int N2 = 0;
    if(fJetsCont2){
        //cout<<"We have a 2nd collection!!"<<endl;
        AliEmcalJet *jet = fJetsCont2->GetNextAcceptJet(0);
        while(jet){
            Float_t NEFpT = jet->Pt()*jet->NEF();
            //cout<<"loop 2 jet 2 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
            N2++;
            Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
            //cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
            fHistJet2->Fill(jetarray);
          //  cout<<" we have a jet! wiht pt = "<<jet->Pt()<<endl;
            jet = fJetsCont2->GetNextAcceptJet();
        }
    }
    //<<"DONE LOOP 2"<<endl;
    
    //cout<<" There are "<<fJetsCont->GetNJets()<<" cont1 jets and "<<fJetsCont2->GetNJets()<<" cont2 jets"<<endl;

    int N1N2 = 0;
    int N1N2m = 0;
   if (fJetsCont&&fJetsCont2) {
         AliEmcalJet *jet1 = fJetsCont->GetNextAcceptJet(0);
         AliEmcalJet *jet2 = fJetsCont2->GetNextAcceptJet(0);
         while(jet1){
             bool ismatched = false;
             Float_t NEFpT1 = jet1->Pt()*jet1->NEF();
             Double_t jetarray1[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),NEFpT1};
             while(jet2){
                 N1N2++;
                 //cout<<"jet 1 has label, pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet1->Pt()<<" "<<jet1->Eta()<<" "<<jet1->Phi()<<" "<<jet1->NEF()<<endl;
                 //<<"jet 2 has lable pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet2->Pt()<<" "<<jet2->Eta()<<" "<<jet2->Phi()<<" "<<jet2->NEF()<<endl;
                 Double_t deta = jet1->Eta()-jet2->Eta();
                 Double_t dphi = RelativePhi(jet1->Phi(),jet2->Phi());
                 Double_t deta2 = deta*deta;
                 Double_t dphi2 = dphi*dphi;
                 Double_t dR = pow(deta2+dphi2,0.5);
                 //cout<<"dR is "<<dR<<endl;
                 Double_t jetarray[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet2->Pt(),jet2->Eta(),jet2->Phi(),jet2->NEF(),dR};
                 if (dR<0.2){
                     N1N2m++;
                     fHistJet1to2->Fill(jetarray);
                     ismatched = true;
                 }
                 jet2 = fJetsCont2->GetNextAcceptJet();
             }
             if (ismatched)
                 fHistJet1m->Fill(jetarray1);
             else
                 fHistJet1nm->Fill(jetarray1);
             jet1 = fJetsCont->GetNextAcceptJet();
             jet2 = fJetsCont2->GetNextAcceptJet(0);
         }
     }

  return kTRUE;
}

//________________________________________________________________________
Float_t AliAnalysisTaskDcalDijetPerf:: RelativePhi(Double_t mphi,Double_t vphi) const
{
    if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
    else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
    if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
    else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
    double dphi = mphi-vphi;
    if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
    else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
    
    return dphi;//dphi in [-Pi, Pi]
}


//________________________________________________________________________
void AliAnalysisTaskDcalDijetPerf::ExecOnce() {

  AliAnalysisTaskEmcalJet::ExecOnce();

  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;

}

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

  return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
}

//________________________________________________________________________
void AliAnalysisTaskDcalDijetPerf::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
b52f1329	1	// $Id$
	2	//
	3	// Dcal dijet performance task
	4	//
	5	// Author: R. Reed
	6
	7	#include <TClonesArray.h>
	8	#include <TH1F.h>
	9	#include <TH2F.h>
	10	#include <TH3F.h>
	11	#include <THnSparse.h>
	12	#include <TList.h>
	13	#include <TLorentzVector.h>
	14
	15	#include "AliVCluster.h"
	16	#include "AliAODCaloCluster.h"
	17	#include "AliESDCaloCluster.h"
	18	#include "AliVTrack.h"
	19	#include "AliEmcalJet.h"
	20	#include "AliRhoParameter.h"
	21	#include "AliLog.h"
	22	#include "AliJetContainer.h"
	23	#include "AliParticleContainer.h"
	24	#include "AliClusterContainer.h"
	25	#include "AliPicoTrack.h"
	26
	27	#include "AliAnalysisTaskDcalDijetPerf.h"
	28
	29	ClassImp(AliAnalysisTaskDcalDijetPerf)
	30
	31	//________________________________________________________________________
	32	AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf() :
	33	AliAnalysisTaskEmcalJet("AliAnalysisTaskDcalDijetPerf", kTRUE),
	34	fHistTracksPt(0),
	35	fHistTracksEtaPhi(0),
	36	fHistClustersPt(0),
	37	fHistLeadingJetPt(0),
	38	fHistJetsPhiEta(0),
	39	fHistJetsPtArea(0),
	40	fHistJetsPtLeadHad(0),
	41	fHistJetsCorrPtArea(0),
	42	fHistJet1(0),
	43	fHistJet1m(0),
	44	fHistJet1nm(0),
	45	fHistJet2(0),
	46	fHistJet1to2(0),
	47	fJetsCont(0),
	48	fJetsCont2(0),
	49	fTracksCont(0),
	50	fCaloClustersCont(0)
	51
	52	{
	53	// Default constructor.
	54
	55	fHistTracksPt = new TH1*[fNcentBins];
	56	fHistTracksEtaPhi = new TH2*[fNcentBins];
	57	fHistClustersPt = new TH1*[fNcentBins];
	58	fHistLeadingJetPt = new TH1*[fNcentBins];
	59	fHistJetsPhiEta = new TH2*[fNcentBins];
	60	fHistJetsPtArea = new TH2*[fNcentBins];
	61	fHistJetsPtLeadHad = new TH2*[fNcentBins];
	62	fHistJetsCorrPtArea = new TH2*[fNcentBins];
	63
	64	for (Int_t i = 0; i < fNcentBins; i++) {
65	fHistTracksPt[i] = 0;
66	fHistTracksEtaPhi[i] = 0;
67	fHistClustersPt[i] = 0;
68	fHistLeadingJetPt[i] = 0;
69	fHistJetsPhiEta[i] = 0;
70	fHistJetsPtArea[i] = 0;
71	fHistJetsPtLeadHad[i] = 0;
72	fHistJetsCorrPtArea[i] = 0;
73	}
74	fHistJet1 = 0;
75	fHistJet1m = 0;
76	fHistJet1nm = 0;
77	fHistJet2 = 0;
78	fHistJet1to2 = 0;
79	SetMakeGeneralHistograms(kTRUE);
80	}
81
82	//________________________________________________________________________
83	AliAnalysisTaskDcalDijetPerf::AliAnalysisTaskDcalDijetPerf(const char *name) :
84	AliAnalysisTaskEmcalJet(name, kTRUE),
85	fHistTracksPt(0),
86	fHistTracksEtaPhi(0),
87	fHistClustersPt(0),
88	fHistLeadingJetPt(0),
89	fHistJetsPhiEta(0),
90	fHistJetsPtArea(0),
91	fHistJetsPtLeadHad(0),
92	fHistJetsCorrPtArea(0),
93	fHistJet1(0),
94	fHistJet1m(0),
95	fHistJet1nm(0),
96	fHistJet2(0),
97	fHistJet1to2(0),
98	fJetsCont(0),
99	fJetsCont2(0),
100	fTracksCont(0),
101	fCaloClustersCont(0)
102	{
103	// Standard constructor.
104
105	fHistTracksPt = new TH1*[fNcentBins];
106	fHistTracksEtaPhi = new TH2*[fNcentBins];
107	fHistClustersPt = new TH1*[fNcentBins];
108	fHistLeadingJetPt = new TH1*[fNcentBins];
109	fHistJetsPhiEta = new TH2*[fNcentBins];
110	fHistJetsPtArea = new TH2*[fNcentBins];
111	fHistJetsPtLeadHad = new TH2*[fNcentBins];
112	fHistJetsCorrPtArea = new TH2*[fNcentBins];
113
114	for (Int_t i = 0; i < fNcentBins; i++) {
115	fHistTracksPt[i] = 0;
116	fHistTracksEtaPhi[i] = 0;
117	fHistClustersPt[i] = 0;
118	fHistLeadingJetPt[i] = 0;
119	fHistJetsPhiEta[i] = 0;
120	fHistJetsPtArea[i] = 0;
121	fHistJetsPtLeadHad[i] = 0;
122	fHistJetsCorrPtArea[i] = 0;
123	}
124	fHistJet1 = 0;
125	fHistJet1m = 0;
126	fHistJet1nm = 0;
127	fHistJet2 = 0;
128	fHistJet1to2 = 0;
129
130	SetMakeGeneralHistograms(kTRUE);
131	}
132
133	//________________________________________________________________________
134	AliAnalysisTaskDcalDijetPerf::~AliAnalysisTaskDcalDijetPerf()
135	{
136	// Destructor.
137	}
138
139	//________________________________________________________________________
140	void AliAnalysisTaskDcalDijetPerf::UserCreateOutputObjects()
141	{
142	// Create user output.
143
144	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
145
146	fJetsCont = GetJetContainer(0);
147	fJetsCont2 = GetJetContainer(1);
148	if(fJetsCont) { //get particles and clusters connected to jets
149	fTracksCont = fJetsCont->GetParticleContainer();
150	fCaloClustersCont = fJetsCont->GetClusterContainer();
151	} else { //no jets, just analysis tracks and clusters
152	fTracksCont = GetParticleContainer(0);
153	fCaloClustersCont = GetClusterContainer(0);
154	}
155	fTracksCont->SetClassName("AliVTrack");
156	//fCaloClustersCont->SetClassName("AliVCluster");
157	TString histname;
158
159	for (Int_t i = 0; i < fNcentBins; i++) {
160	if (fParticleCollArray.GetEntriesFast()>0) {
161	histname = "fHistTracksPt_";
162	histname += i;
163	fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
164	fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
165	fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
166	fOutput->Add(fHistTracksPt[i]);
167	histname = "fHistTracksEtaPhi_";
168	histname += i;
169	fHistTracksEtaPhi[i] = new TH2F(histname.Data(), histname.Data(), fNbins, -0.7, 0.7, fNbins, 0, TMath::Pi()*2);
170	fHistTracksEtaPhi[i]->GetXaxis()->SetTitle("eta");
171	fHistTracksEtaPhi[i]->GetYaxis()->SetTitle("phi");
172	fOutput->Add(fHistTracksEtaPhi[i]);
173	}
174
175	if (fClusterCollArray.GetEntriesFast()>0) {
176	histname = "fHistClustersPt_";
177	histname += i;
178	fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
179	fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
180	fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
181	fOutput->Add(fHistClustersPt[i]);
182	}
183
184	if (fJetCollArray.GetEntriesFast()>0) {
185	histname = "fHistLeadingJetPt_";
186	histname += i;
187	fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
188	fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
189	fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
190	fOutput->Add(fHistLeadingJetPt[i]);
191
192	histname = "fHistJetsPhiEta_";
193	histname += i;
194	fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
195	fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
196	fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
197	fOutput->Add(fHistJetsPhiEta[i]);
198
199	histname = "fHistJetsPtArea_";
200	histname += i;
201	fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
202	fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
203	fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
204	fOutput->Add(fHistJetsPtArea[i]);
205
206	histname = "fHistJetsPtLeadHad_";
207	histname += i;
208	fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
209	fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
210	fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
211	fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
212	fOutput->Add(fHistJetsPtLeadHad[i]);
213
214	if (!(GetJetContainer()->GetRhoName().IsNull())) {
215	histname = "fHistJetsCorrPtArea_";
216	histname += i;
217	fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
218	fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
219	fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
220	fOutput->Add(fHistJetsCorrPtArea[i]);
221	}
222	}
223	}
224
225	Int_t nbins[] = {150,100,100,100,150};
226	Double_t xmin[] = { 0,-0.7, 0, 0, 0};
227	Double_t xmax[] = {150, 0.7,TMath::TwoPi(), 1, 150};
228	fHistJet1 = new THnSparseF("Jets1Collection","Jets1Collection",5,nbins,xmin,xmax);
229	fOutput->Add(fHistJet1);
230	fHistJet1m = new THnSparseF("Jets1CollectionMatched","Jets1Collection",5,nbins,xmin,xmax);
231	fOutput->Add(fHistJet1m);
232	fHistJet1nm = new THnSparseF("Jets1CollectionNotMatched","Jets1Collection",5,nbins,xmin,xmax);
233	fOutput->Add(fHistJet1nm);
234	fHistJet2 = new THnSparseF("Jets2Collection","Jets2Collection",5,nbins,xmin,xmax);
235	fOutput->Add(fHistJet2);
236
237	Int_t nbins2[] = {150,100,100,100,150,100,100,100,100};
238	Double_t xmin2[] = {0,-0.7,0,0,0,-0.7,0,0,0};
239	Double_t xmax2[] = {150,0.7,6.28,1,150,0.7,6.28,1,0.2};
240	fHistJet1to2 = new THnSparseF("Jets1to2Collection","Jets1to2Collection",9,nbins2,xmin2,xmax2);
241	fOutput->Add(fHistJet1to2);
242
243	PostData(1, fOutput); // Post data for ALL output slots > 0 here.
244	}
245
246	//________________________________________________________________________
247	Bool_t AliAnalysisTaskDcalDijetPerf::FillHistograms()
248	{
249	// Fill histograms.
250
251	if (fTracksCont) {
252	AliVTrack track = static_cast<AliVTrack>(fTracksCont->GetNextAcceptParticle(0));
253	while(track) {
254	fHistTracksPt[fCentBin]->Fill(track->Pt());
255	fHistTracksEtaPhi[fCentBin]->Fill(track->Eta(),track->Phi());
256	track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
257	}
258	}
259
260	if (fCaloClustersCont) {
261	AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
262	while(cluster) {
263	TLorentzVector nPart;
264	cluster->GetMomentum(nPart, fVertex);
265	fHistClustersPt[fCentBin]->Fill(nPart.Pt());
266
267	cluster = fCaloClustersCont->GetNextAcceptCluster();
268	}
269	}
270
271	int N1 = 0;
272	if (fJetsCont) {
273	AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
274	while(jet) {
275	Float_t NEFpT = jet->Pt()*jet->NEF();
276	N1++;
277	//cout<<"loop 1 jet 1 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<" and NEFpT "<<NEFpT<<endl;
278	fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
279	fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
280	Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
281	fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
282
283	//110 degrees in the azimuthal angle \|eta\|<0.7
284	//EMCal Tower 0.0143 x 0.0143
285	//16x16 = 0.2288 x 0.2288 -> R = 0.1144
286
287	Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
288	//cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
289	fHistJet1->Fill(jetarray);
290	if (fHistJetsCorrPtArea[fCentBin]) {
291	Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
292	fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
293	}
294	jet = fJetsCont->GetNextAcceptJet();
295	}
296
297	jet = fJetsCont->GetLeadingJet();
298	if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
299	}
300	// cout<<"DONE LOOP 1"<<endl;
301	int N2 = 0;
302	if(fJetsCont2){
303	//cout<<"We have a 2nd collection!!"<<endl;
304	AliEmcalJet *jet = fJetsCont2->GetNextAcceptJet(0);
305	while(jet){
306	Float_t NEFpT = jet->Pt()*jet->NEF();
307	//cout<<"loop 2 jet 2 has label, pt,eta,phi,NEF = "<<jet->GetLabel()<<" "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
308	N2++;
309	Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
310	//cout<<"Filling "<<jet->Pt()<<" "<<jet->Eta()<<" "<<jet->Phi()<<" "<<jet->NEF()<<endl;
311	fHistJet2->Fill(jetarray);
312	// cout<<" we have a jet! wiht pt = "<<jet->Pt()<<endl;
313	jet = fJetsCont2->GetNextAcceptJet();
314	}
315	}
316	//<<"DONE LOOP 2"<<endl;
317
318	//cout<<" There are "<<fJetsCont->GetNJets()<<" cont1 jets and "<<fJetsCont2->GetNJets()<<" cont2 jets"<<endl;
319
320	int N1N2 = 0;
321	int N1N2m = 0;
322	if (fJetsCont&&fJetsCont2) {
323	AliEmcalJet *jet1 = fJetsCont->GetNextAcceptJet(0);
324	AliEmcalJet *jet2 = fJetsCont2->GetNextAcceptJet(0);
325	while(jet1){
326	bool ismatched = false;
327	Float_t NEFpT1 = jet1->Pt()*jet1->NEF();
328	Double_t jetarray1[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),NEFpT1};
329	while(jet2){
330	N1N2++;
331	//cout<<"jet 1 has label, pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet1->Pt()<<" "<<jet1->Eta()<<" "<<jet1->Phi()<<" "<<jet1->NEF()<<endl;
332	//<<"jet 2 has lable pt,eta,phi,NEF = "<<jet1->GetLabel()<<" "<<jet2->Pt()<<" "<<jet2->Eta()<<" "<<jet2->Phi()<<" "<<jet2->NEF()<<endl;
333	Double_t deta = jet1->Eta()-jet2->Eta();
334	Double_t dphi = RelativePhi(jet1->Phi(),jet2->Phi());
335	Double_t deta2 = deta*deta;
336	Double_t dphi2 = dphi*dphi;
337	Double_t dR = pow(deta2+dphi2,0.5);
338	//cout<<"dR is "<<dR<<endl;
339	Double_t jetarray[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet2->Pt(),jet2->Eta(),jet2->Phi(),jet2->NEF(),dR};
340	if (dR<0.2){
341	N1N2m++;
342	fHistJet1to2->Fill(jetarray);
343	ismatched = true;
344	}
345	jet2 = fJetsCont2->GetNextAcceptJet();
346	}
347	if (ismatched)
348	fHistJet1m->Fill(jetarray1);
349	else
350	fHistJet1nm->Fill(jetarray1);
351	jet1 = fJetsCont->GetNextAcceptJet();
352	jet2 = fJetsCont2->GetNextAcceptJet(0);
353	}
354	}
355
356	return kTRUE;
357	}
358
359	//________________________________________________________________________
360	Float_t AliAnalysisTaskDcalDijetPerf:: RelativePhi(Double_t mphi,Double_t vphi) const
361	{
362	if (vphi < -1TMath::Pi()) vphi += (2TMath::Pi());
363	else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
364	if (mphi < -1TMath::Pi()) mphi += (2TMath::Pi());
365	else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
366	double dphi = mphi-vphi;
367	if (dphi < -1TMath::Pi()) dphi += (2TMath::Pi());
368	else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
369
370	return dphi;//dphi in [-Pi, Pi]
371	}
372
373
374	//________________________________________________________________________
375	void AliAnalysisTaskDcalDijetPerf::ExecOnce() {
376
377	AliAnalysisTaskEmcalJet::ExecOnce();
378
379	if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
380	if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
381	if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
382
383	}
384
385	//________________________________________________________________________
386	Bool_t AliAnalysisTaskDcalDijetPerf::Run()
387	{
388	// Run analysis code here, if needed. It will be executed before FillHistograms().
389
390	return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
391	}
392
393	//________________________________________________________________________
394	void AliAnalysisTaskDcalDijetPerf::Terminate(Option_t *)
395	{
396	// Called once at the end of the analysis.
397	}