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

#include "AliAnalysisTaskEmcalJetHadCorQA.h"

#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THnSparse.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TParameter.h>
#include <TParticle.h>
#include <TTree.h>
#include <TVector3.h>

#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTask.h"
#include "AliCentrality.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliEmcalJet.h"
#include "AliVCluster.h"
#include "AliVTrack.h"
#include "AliRhoParameter.h"
#include "AliEmcalParticle.h"
#include "AliPicoTrack.h"
#include "AliEMCALGeometry.h"
using std::vector;

ClassImp(AliAnalysisTaskEmcalJetHadCorQA)

//________________________________________________________________________
AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() : 
  AliAnalysisTaskEmcalJet("jethadcor",kFALSE), 
  fCalo2Name(),
  fCaloClusters2(),
  fMCParticlesName(),
  fMCParticles(),
  fHistRhovsCent(0),
  fHistNjetvsCent(0),
  fHistNTMatchvsPtvsNtack0(0)
{
  for (int i = 0;i<3;i++){
  fHistNEFvsPt[i] = 0;
  fHistNTMatchvsPt[i] = 0;  
  fHistNCMatchvsPt[i] = 0;
  fHistHadCorvsPt[i] = 0;
  fHistNEFvsPtBias[i] = 0;
  fHistNconvsPt[i] = 0;  
  fHistNtvsPt[i] = 0;    
  fHistNcvsPt[i] = 0;    
  fHistNTMatchvsPtBias[i] = 0;  
  fHistNCMatchvsPtBias[i] = 0;
  fHistHadCorvsPtBias[i] = 0;
  fHistNconvsPtBias[i] = 0;  
  fHistNtvsPtBias[i] = 0;    
  fHistNcvsPtBias[i] = 0;    
  }
  // Default constructor.
 
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
  AliAnalysisTaskEmcalJet(name,kTRUE),
  fCalo2Name(),
  fCaloClusters2(),
  fMCParticlesName(),
  fMCParticles(),
  fHistRhovsCent(0),
  fHistNjetvsCent(0),
  fHistNTMatchvsPtvsNtack0(0)
 { 
  for (int i = 0;i<3;i++){
  fHistNEFvsPt[i] = 0;
  fHistNTMatchvsPt[i] = 0;  
  fHistNCMatchvsPt[i] = 0;
  fHistHadCorvsPt[i] = 0;
  fHistNEFvsPtBias[i] = 0;
  fHistNconvsPt[i] = 0;  
  fHistNtvsPt[i] = 0;    
  fHistNcvsPt[i] = 0;    
  fHistNTMatchvsPtBias[i] = 0;  
  fHistNCMatchvsPtBias[i] = 0;
  fHistHadCorvsPtBias[i] = 0;
  fHistNconvsPtBias[i] = 0;  
  fHistNtvsPtBias[i] = 0;    
  fHistNcvsPtBias[i] = 0;    
  }
   SetMakeGeneralHistograms(kTRUE);
 }

//________________________________________________________________________
void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
{
  if (! fCreateHisto)
    return;
  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  fHistRhovsCent             = new TH2F("RhovsCent",              "RhovsCent",             100, 0.0, 100.0, 500, 0, 500);
  fHistNjetvsCent            = new TH2F("NjetvsCent",             "NjetvsCent",            100, 0.0, 100.0, 100, 0, 100);

  for (int i = 0;i<3;i++){
    char name[200];
    TString title;
    sprintf(name,"NEFvsPt%i",i);
    fHistNEFvsPt[i]         = new TH2F(name, name, 100,0,1,500,-250,250);
    fOutput->Add(fHistNEFvsPt[i]);
    sprintf(name,"NTMatchvsPt%i",i);
    fHistNTMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);  
    fOutput->Add(fHistNTMatchvsPt[i]);
    sprintf(name,"NCMatchvsPt%i",i);
    fHistNCMatchvsPt[i]     = new TH2F(name, name, 100,0,100,500,-250,250);
    fOutput->Add(fHistNCMatchvsPt[i]);
    sprintf(name,"HadCorvsPt%i",i);
    fHistHadCorvsPt[i]      = new TH2F(name, name, 1000,0,500,500,-250,250);
    fOutput->Add(fHistHadCorvsPt[i]);
    sprintf(name,"NconvsPt%i",i);
    fHistNconvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNconvsPt[i]);
    sprintf(name,"NtvsPt%i",i);
    fHistNtvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNtvsPt[i]);
    sprintf(name,"NcvsPt%i",i);
    fHistNcvsPt[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNcvsPt[i]);
    sprintf(name,"NEFvsPtBias%i",i);
    fHistNEFvsPtBias[i]     = new TH2F(name, name, 100,0,1,500,-250,250);
    fOutput->Add(fHistNEFvsPtBias[i]);
    sprintf(name,"NTMatchvsPtBias%i",i);
    fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);  
    fOutput->Add(fHistNTMatchvsPtBias[i]);
    sprintf(name,"NCMatchvsPtBias%i",i);
    fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
    fOutput->Add(fHistNCMatchvsPtBias[i]);
    sprintf(name,"HadCorvsPtBias%i",i);
    fHistHadCorvsPtBias[i]  = new TH2F(name, name, 1000,0,500,500,-250,250);
    fOutput->Add(fHistHadCorvsPtBias[i]);
   sprintf(name,"NconvsPtBias%i",i);
    fHistNconvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNconvsPtBias[i]);
    sprintf(name,"NtvsPtBias%i",i);
    fHistNtvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNtvsPtBias[i]);
    sprintf(name,"NcvsPtBias%i",i);
    fHistNcvsPtBias[i]      = new TH2F(name, name, 200,0,200,500,-250,250);
    fOutput->Add(fHistNcvsPtBias[i]);
  }
  fHistNTMatchvsPtvsNtack0   = new TH3F("NTMmatchvsPtvsNtrack0",  "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);

  // for (Int_t i = 0;i<6;++i){
  //   name = TString(Form("JetPtvsTrackPt_%i",i));
  //   title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
  //   fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
  //   fOutput->Add(fHistJetPtvsTrackPt[i]);
   
  // }
  fOutput->Add(fHistRhovsCent);
  fOutput->Add(fHistNjetvsCent);
  fOutput->Add(fHistNTMatchvsPtvsNtack0);
   PostData(1, fOutput);
}

//________________________________________________________________________

Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const 
{
  // Get centrality bin.

  Int_t centbin = -1;
  if (cent>=0 && cent<10)
    centbin = 0;
  else if (cent>=10 && cent<30)
    centbin = 1;
  else if (cent>=30 && cent<50)
    centbin = 2;
  else if (cent>50)
    centbin =3;
  return centbin;
}

//________________________________________________________________________

Float_t AliAnalysisTaskEmcalJetHadCorQA:: 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 AliAnalysisTaskEmcalJetHadCorQA::ExecOnce()
{
  AliAnalysisTaskEmcalJet::ExecOnce();

//   AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
//   if (fCaloName == "CaloClusters")
//     am->LoadBranch("CaloClusters");

  if (!fCalo2Name.IsNull() && !fCaloClusters2){
    fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
    if (!fCaloClusters2){
      AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
      fInitialized = kFALSE;
      return;
    }
  }

  if (!fMCParticlesName.IsNull() && !fMCParticles){
    fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
    if (!fMCParticles){
      AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
      fInitialized = kFALSE;
      return;
    }
  }
    //TString             fMCParticlesName;
    // TClonesArray       *fMCParticles;
    //  fCaloClusters2(),
 //    else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
//       AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
//       fCaloClusters2 = 0;
//       fInitialized = kFALSE;
//       return;
    //  }
}


//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
{  
  Int_t centbin = GetCentBin(fCent);
  //for pp analyses we will just use the first centrality bin
  if (GetBeamType()==0)
    centbin = 0;
  if (centbin>2)
    return kTRUE;
  if (!fTracks)
    return kTRUE;
  if (!fCaloClusters)
    return kTRUE;
  if (!fCaloClusters2)
    return kTRUE;
  const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
  const Int_t nTrack = fTracks->GetEntriesFast();
  //  const Int_t nMC = fMCParticles->GetEntriesFast();

  TString fRhoScaledName = fRhoName;
  fRho = GetRhoFromEvent(fRhoScaledName);
  fRhoVal = fRho->GetVal();
  if (GetBeamType()==0)
    fRhoVal = 0;
  fHistRhovsCent->Fill(fCent,fRhoVal);
  const Int_t Njets = fJets->GetEntriesFast();

  Int_t NjetAcc = 0;

 
  for (Int_t iJets = 0; iJets < Njets; ++iJets) {
    Int_t TrackMatch = 0;
     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
     if (!jet)
       continue; 
     if (jet->Area()==0)
       continue;
     if (jet->Pt()<0.1)
       continue;
     if (jet->MaxTrackPt()>100)
       continue;
     if (! AcceptJet(jet))
       continue;
     NjetAcc++;
     vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
     vector<Int_t> cluster_id2;
     Double_t Esub = 0; //total E subtracted from the jet
     Double_t jetPt = -500;
     jetPt = jet->Pt()-jet->Area()*fRhoVal;    
     for (int i = 0;i<nTrack;i++){
       AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(i));
       if (!emctrack)
	 continue;
       if (!emctrack->IsEMCAL())
	 continue;
       AliVTrack *track = (AliVTrack*)emctrack->GetTrack();
       if (! track)
	 continue;
       if (! AcceptTrack(track))
	 continue;
       if (! IsJetTrack(jet,i,false))
	 continue;
       Int_t iClus = track->GetEMCALcluster();
       if (iClus<0)
	 continue;
       //we have the id of the matched cluster of a track constituent
       //       cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
       bool ischecked = false;
       for (UInt_t icid = 0;icid<cluster_id.size();icid++)
	 if (cluster_id[icid] == iClus)
	   ischecked = true; // we've already looked at this uncorrected cluster
       if (ischecked)
	 continue; //no need to go further
       AliEmcalParticle *emcluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
       AliVCluster *cluster = emcluster->GetCluster();
       if (! cluster)
	 continue;
       if (! AcceptCluster(cluster))
	 continue;
       Double_t etadiff = 999;
       Double_t phidiff = 999;
       AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
       if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
	 continue;
       TrackMatch++; //this cluster has been matched, let's add it to the list
       cluster_id.push_back(iClus);
       Int_t ismatch = -1;
       //now we need to find its matched corrected cluster if any
       for (Int_t ic = 0;ic < nCluster2;ic++){
	 //we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
	 AliVCluster *emcluster2 = static_cast<AliVCluster*>(fCaloClusters2->At(ic));
	 if (! emcluster2)
	   continue;
	 if (! AcceptCluster(emcluster2))
	   continue;
	 if (! IsJetCluster(jet,ic,false))
	   continue;
	 TLorentzVector nPart;
	 cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
	 TLorentzVector nPart2;
	 emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
	 float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
	 if (R < 0.001){//this cluster stayed in the jet!
	   cluster_id2.push_back(ic);
	   ismatch = ic;
	 }
       } // end of cluster2 loop
       //we only get here if the track was matched to a cluster that hadn't been looked at before
       if (ismatch < 0) //this cluster was entirely deleted
	 Esub+=cluster->E();
       else{ //get the corrected cluster
	 AliVCluster *emcluster2temp = static_cast<AliVCluster*>(fCaloClusters2->At(ismatch));
	 Esub+=(cluster->E() - emcluster2temp->E());
       }
       
     } // end of track loop
     fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
     fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
     fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
     fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
     fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
     fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);   
     fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);   
     if (jet->MaxTrackPt()<5.0)
       continue;
     fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
     fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
     fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
     fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
     fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);   
     fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);  
     fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt); 
     if (centbin == 0)
       fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
  }
  fHistNjetvsCent->Fill(fCent,NjetAcc);
  return kTRUE;
}
Commit	Line	Data
6d9675ae	1	#include "AliAnalysisTaskEmcalJetHadCorQA.h"
	2
	3	#include <TCanvas.h>
	4	#include <TChain.h>
	5	#include <TClonesArray.h>
	6	#include <TH1F.h>
	7	#include <TH2F.h>
	8	#include <TH3F.h>
	9	#include <THnSparse.h>
	10	#include <TList.h>
	11	#include <TLorentzVector.h>
	12	#include <TParameter.h>
	13	#include <TParticle.h>
	14	#include <TTree.h>
	15	#include <TVector3.h>
	16
	17	#include "AliAODEvent.h"
	18	#include "AliAnalysisManager.h"
	19	#include "AliAnalysisTask.h"
	20	#include "AliCentrality.h"
	21	#include "AliESDEvent.h"
	22	#include "AliESDInputHandler.h"
	23	#include "AliEmcalJet.h"
	24	#include "AliVCluster.h"
	25	#include "AliVTrack.h"
	26	#include "AliRhoParameter.h"
	27	#include "AliEmcalParticle.h"
	28	#include "AliPicoTrack.h"
	29	#include "AliEMCALGeometry.h"
c2aad3ae	30	using std::vector;
6d9675ae	31
	32	ClassImp(AliAnalysisTaskEmcalJetHadCorQA)
	33
	34	//________________________________________________________________________
	35	AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() :
9239b066	36	AliAnalysisTaskEmcalJet("jethadcor",kFALSE),
6d9675ae	37	fCalo2Name(),
6d9675ae	38	fCaloClusters2(),
ba03cebd	39	fMCParticlesName(),
ba03cebd	40	fMCParticles(),
6d9675ae	41	fHistRhovsCent(0),
fb99636c	42	fHistNjetvsCent(0),
fb99636c	43	fHistNTMatchvsPtvsNtack0(0)
6d9675ae	44	{
	45	for (int i = 0;i<3;i++){
	46	fHistNEFvsPt[i] = 0;
	47	fHistNTMatchvsPt[i] = 0;
	48	fHistNCMatchvsPt[i] = 0;
	49	fHistHadCorvsPt[i] = 0;
	50	fHistNEFvsPtBias[i] = 0;
ba03cebd	51	fHistNconvsPt[i] = 0;
	52	fHistNtvsPt[i] = 0;
	53	fHistNcvsPt[i] = 0;
6d9675ae	54	fHistNTMatchvsPtBias[i] = 0;
	55	fHistNCMatchvsPtBias[i] = 0;
	56	fHistHadCorvsPtBias[i] = 0;
ba03cebd	57	fHistNconvsPtBias[i] = 0;
	58	fHistNtvsPtBias[i] = 0;
	59	fHistNcvsPtBias[i] = 0;
6d9675ae	60	}
	61	// Default constructor.
	62
	63	SetMakeGeneralHistograms(kTRUE);
	64	}
	65
	66	//________________________________________________________________________
	67	AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
9239b066	68	AliAnalysisTaskEmcalJet(name,kTRUE),
6d9675ae	69	fCalo2Name(),
6d9675ae	70	fCaloClusters2(),
ba03cebd	71	fMCParticlesName(),
ba03cebd	72	fMCParticles(),
6d9675ae	73	fHistRhovsCent(0),
fb99636c	74	fHistNjetvsCent(0),
fb99636c	75	fHistNTMatchvsPtvsNtack0(0)
6d9675ae	76	{
	77	for (int i = 0;i<3;i++){
	78	fHistNEFvsPt[i] = 0;
	79	fHistNTMatchvsPt[i] = 0;
	80	fHistNCMatchvsPt[i] = 0;
	81	fHistHadCorvsPt[i] = 0;
	82	fHistNEFvsPtBias[i] = 0;
ba03cebd	83	fHistNconvsPt[i] = 0;
	84	fHistNtvsPt[i] = 0;
	85	fHistNcvsPt[i] = 0;
6d9675ae	86	fHistNTMatchvsPtBias[i] = 0;
	87	fHistNCMatchvsPtBias[i] = 0;
	88	fHistHadCorvsPtBias[i] = 0;
ba03cebd	89	fHistNconvsPtBias[i] = 0;
	90	fHistNtvsPtBias[i] = 0;
	91	fHistNcvsPtBias[i] = 0;
6d9675ae	92	}
	93	SetMakeGeneralHistograms(kTRUE);
	94	}
	95
	96	//________________________________________________________________________
	97	void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
	98	{
	99	if (! fCreateHisto)
	100	return;
9239b066	101	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
6d9675ae	102	fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 500, 0, 500);
	103	fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
	104
	105	for (int i = 0;i<3;i++){
	106	char name[200];
	107	TString title;
	108	sprintf(name,"NEFvsPt%i",i);
	109	fHistNEFvsPt[i] = new TH2F(name, name, 100,0,1,500,-250,250);
	110	fOutput->Add(fHistNEFvsPt[i]);
	111	sprintf(name,"NTMatchvsPt%i",i);
	112	fHistNTMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
	113	fOutput->Add(fHistNTMatchvsPt[i]);
	114	sprintf(name,"NCMatchvsPt%i",i);
	115	fHistNCMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
	116	fOutput->Add(fHistNCMatchvsPt[i]);
	117	sprintf(name,"HadCorvsPt%i",i);
	118	fHistHadCorvsPt[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
	119	fOutput->Add(fHistHadCorvsPt[i]);
ba03cebd	120	sprintf(name,"NconvsPt%i",i);
	121	fHistNconvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	122	fOutput->Add(fHistNconvsPt[i]);
	123	sprintf(name,"NtvsPt%i",i);
	124	fHistNtvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	125	fOutput->Add(fHistNtvsPt[i]);
	126	sprintf(name,"NcvsPt%i",i);
	127	fHistNcvsPt[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	128	fOutput->Add(fHistNcvsPt[i]);
6d9675ae	129	sprintf(name,"NEFvsPtBias%i",i);
	130	fHistNEFvsPtBias[i] = new TH2F(name, name, 100,0,1,500,-250,250);
	131	fOutput->Add(fHistNEFvsPtBias[i]);
	132	sprintf(name,"NTMatchvsPtBias%i",i);
	133	fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
	134	fOutput->Add(fHistNTMatchvsPtBias[i]);
	135	sprintf(name,"NCMatchvsPtBias%i",i);
	136	fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
	137	fOutput->Add(fHistNCMatchvsPtBias[i]);
	138	sprintf(name,"HadCorvsPtBias%i",i);
	139	fHistHadCorvsPtBias[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
	140	fOutput->Add(fHistHadCorvsPtBias[i]);
ba03cebd	141	sprintf(name,"NconvsPtBias%i",i);
	142	fHistNconvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	143	fOutput->Add(fHistNconvsPtBias[i]);
	144	sprintf(name,"NtvsPtBias%i",i);
	145	fHistNtvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	146	fOutput->Add(fHistNtvsPtBias[i]);
	147	sprintf(name,"NcvsPtBias%i",i);
	148	fHistNcvsPtBias[i] = new TH2F(name, name, 200,0,200,500,-250,250);
	149	fOutput->Add(fHistNcvsPtBias[i]);
6d9675ae	150	}
	151	fHistNTMatchvsPtvsNtack0 = new TH3F("NTMmatchvsPtvsNtrack0", "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
	152
	153	// for (Int_t i = 0;i<6;++i){
	154	// name = TString(Form("JetPtvsTrackPt_%i",i));
	155	// title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
	156	// fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
	157	// fOutput->Add(fHistJetPtvsTrackPt[i]);
	158
	159	// }
	160	fOutput->Add(fHistRhovsCent);
	161	fOutput->Add(fHistNjetvsCent);
	162	fOutput->Add(fHistNTMatchvsPtvsNtack0);
	163	PostData(1, fOutput);
	164	}
	165
	166	//________________________________________________________________________
	167
	168	Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const
	169	{
	170	// Get centrality bin.
	171
	172	Int_t centbin = -1;
	173	if (cent>=0 && cent<10)
	174	centbin = 0;
	175	else if (cent>=10 && cent<30)
	176	centbin = 1;
	177	else if (cent>=30 && cent<50)
	178	centbin = 2;
	179	else if (cent>50)
	180	centbin =3;
	181	return centbin;
	182	}
	183
	184	//________________________________________________________________________
	185
	186	Float_t AliAnalysisTaskEmcalJetHadCorQA:: RelativePhi(Double_t mphi,Double_t vphi) const
	187	{
	188	if (vphi < -1TMath::Pi()) vphi += (2TMath::Pi());
	189	else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
	190	if (mphi < -1TMath::Pi()) mphi += (2TMath::Pi());
	191	else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
	192	double dphi = mphi-vphi;
	193	if (dphi < -1TMath::Pi()) dphi += (2TMath::Pi());
	194	else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
	195
	196	return dphi;//dphi in [-Pi, Pi]
	197	}
	198
	199	//________________________________________________________________________
9239b066	200	void AliAnalysisTaskEmcalJetHadCorQA::ExecOnce()
	201	{
	202	AliAnalysisTaskEmcalJet::ExecOnce();
6d9675ae	203
	204	// AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
	205	// if (fCaloName == "CaloClusters")
	206	// am->LoadBranch("CaloClusters");
	207
	208	if (!fCalo2Name.IsNull() && !fCaloClusters2){
	209	fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
	210	if (!fCaloClusters2){
	211	AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
	212	fInitialized = kFALSE;
	213	return;
	214	}
ba03cebd	215	}
	216
	217	if (!fMCParticlesName.IsNull() && !fMCParticles){
	218	fMCParticles = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticlesName));
	219	if (!fMCParticles){
	220	AliError(Form("%s: Could not retrieve MC Particles %s!",GetName(),fMCParticlesName.Data()));
	221	fInitialized = kFALSE;
	222	return;
	223	}
	224	}
	225	//TString fMCParticlesName;
	226	// TClonesArray *fMCParticles;
6d9675ae	227	// fCaloClusters2(),
	228	// else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
	229	// AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
	230	// fCaloClusters2 = 0;
	231	// fInitialized = kFALSE;
	232	// return;
	233	// }
6d9675ae	234	}
	235
	236
	237	//________________________________________________________________________
	238	Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
	239	{
	240	Int_t centbin = GetCentBin(fCent);
	241	//for pp analyses we will just use the first centrality bin
ba03cebd	242	if (GetBeamType()==0)
6d9675ae	243	centbin = 0;
	244	if (centbin>2)
	245	return kTRUE;
	246	if (!fTracks)
	247	return kTRUE;
ba03cebd	248	if (!fCaloClusters)
	249	return kTRUE;
	250	if (!fCaloClusters2)
	251	return kTRUE;
6d9675ae	252	const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
6d9675ae	253	const Int_t nTrack = fTracks->GetEntriesFast();
ba03cebd	254	// const Int_t nMC = fMCParticles->GetEntriesFast();
6d9675ae	255
	256	TString fRhoScaledName = fRhoName;
	257	fRho = GetRhoFromEvent(fRhoScaledName);
	258	fRhoVal = fRho->GetVal();
ba03cebd	259	if (GetBeamType()==0)
ba03cebd	260	fRhoVal = 0;
6d9675ae	261	fHistRhovsCent->Fill(fCent,fRhoVal);
	262	const Int_t Njets = fJets->GetEntriesFast();
	263
	264	Int_t NjetAcc = 0;
	265
	266
	267	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
	268	Int_t TrackMatch = 0;
	269	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(iJets));
	270	if (!jet)
	271	continue;
	272	if (jet->Area()==0)
	273	continue;
	274	if (jet->Pt()<0.1)
	275	continue;
	276	if (jet->MaxTrackPt()>100)
	277	continue;
	278	if (! AcceptJet(jet))
	279	continue;
	280	NjetAcc++;
	281	vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
	282	vector<Int_t> cluster_id2;
	283	Double_t Esub = 0; //total E subtracted from the jet
	284	Double_t jetPt = -500;
	285	jetPt = jet->Pt()-jet->Area()*fRhoVal;
	286	for (int i = 0;i<nTrack;i++){
	287	AliEmcalParticle emctrack = static_cast<AliEmcalParticle>(fTracks->At(i));
	288	if (!emctrack)
	289	continue;
	290	if (!emctrack->IsEMCAL())
	291	continue;
	292	AliVTrack track = (AliVTrack)emctrack->GetTrack();
	293	if (! track)
	294	continue;
	295	if (! AcceptTrack(track))
	296	continue;
	297	if (! IsJetTrack(jet,i,false))
	298	continue;
	299	Int_t iClus = track->GetEMCALcluster();
	300	if (iClus<0)
	301	continue;
ba03cebd	302	//we have the id of the matched cluster of a track constituent
ba03cebd	303	// cout<<"track label for matched,accepted jet track is "<<track->GetLabel()<<endl;
6d9675ae	304	bool ischecked = false;
fb99636c	305	for (UInt_t icid = 0;icid<cluster_id.size();icid++)
6d9675ae	306	if (cluster_id[icid] == iClus)
	307	ischecked = true; // we've already looked at this uncorrected cluster
	308	if (ischecked)
	309	continue; //no need to go further
	310	AliEmcalParticle emcluster = static_cast<AliEmcalParticle>(fCaloClusters->At(iClus));
	311	AliVCluster *cluster = emcluster->GetCluster();
	312	if (! cluster)
	313	continue;
	314	if (! AcceptCluster(cluster))
	315	continue;
	316	Double_t etadiff = 999;
	317	Double_t phidiff = 999;
	318	AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
	319	if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
	320	continue;
	321	TrackMatch++; //this cluster has been matched, let's add it to the list
	322	cluster_id.push_back(iClus);
	323	Int_t ismatch = -1;
	324	//now we need to find its matched corrected cluster if any
	325	for (Int_t ic = 0;ic < nCluster2;ic++){
	326	//we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
	327	AliVCluster emcluster2 = static_cast<AliVCluster>(fCaloClusters2->At(ic));
	328	if (! emcluster2)
	329	continue;
	330	if (! AcceptCluster(emcluster2))
	331	continue;
	332	if (! IsJetCluster(jet,ic,false))
	333	continue;
	334	TLorentzVector nPart;
	335	cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
	336	TLorentzVector nPart2;
	337	emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
	338	float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
	339	if (R < 0.001){//this cluster stayed in the jet!
	340	cluster_id2.push_back(ic);
	341	ismatch = ic;
	342	}
	343	} // end of cluster2 loop
	344	//we only get here if the track was matched to a cluster that hadn't been looked at before
	345	if (ismatch < 0) //this cluster was entirely deleted
	346	Esub+=cluster->E();
	347	else{ //get the corrected cluster
	348	AliVCluster emcluster2temp = static_cast<AliVCluster>(fCaloClusters2->At(ismatch));
	349	Esub+=(cluster->E() - emcluster2temp->E());
	350	}
	351
	352	} // end of track loop
	353	fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
	354	fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
ba03cebd	355	fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
6d9675ae	356	fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
ba03cebd	357	fHistNconvsPt[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
	358	fHistNtvsPt[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
	359	fHistNcvsPt[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
6d9675ae	360	if (jet->MaxTrackPt()<5.0)
	361	continue;
	362	fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
	363	fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
	364	fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
ba03cebd	365	fHistNconvsPtBias[centbin]->Fill(jet->GetNumberOfConstituents(),jetPt);
	366	fHistNtvsPtBias[centbin]->Fill(jet->GetNumberOfTracks(),jetPt);
	367	fHistNcvsPtBias[centbin]->Fill(jet->GetNumberOfClusters(),jetPt);
	368	fHistNCMatchvsPt[centbin]->Fill(cluster_id.size(),jetPt);
6d9675ae	369	if (centbin == 0)
ba03cebd	370	fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,jet->GetNumberOfTracks());
6d9675ae	371	}
	372	fHistNjetvsCent->Fill(fCent,NjetAcc);
	373	return kTRUE;
	374	}
	375
	376
	377
	378
	379