[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"

ClassImp(AliAnalysisTaskEmcalJetHadCorQA)

//________________________________________________________________________
AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() : 
  AliAnalysisTaskEmcalJet("spectra",kFALSE), 
  fCalo2Name(),
  fCaloClusters2(),
  fHistRhovsCent(0),
  fHistNjetvsCent(0)
{
  for (int i = 0;i<3;i++){
  fHistNEFvsPt[i] = 0;
  fHistNTMatchvsPt[i] = 0;  
  fHistNCMatchvsPt[i] = 0;
  fHistHadCorvsPt[i] = 0;
  fHistNEFvsPtBias[i] = 0;
  fHistNTMatchvsPtBias[i] = 0;  
  fHistNCMatchvsPtBias[i] = 0;
  fHistHadCorvsPtBias[i] = 0;
  }
  // Default constructor.
 
  SetMakeGeneralHistograms(kTRUE);
}

//________________________________________________________________________
AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
  AliAnalysisTaskEmcalJet(name,kTRUE),
  fCalo2Name(),
  fCaloClusters2(),
  fHistRhovsCent(0),
  fHistNjetvsCent(0)
 { 
  for (int i = 0;i<3;i++){
  fHistNEFvsPt[i] = 0;
  fHistNTMatchvsPt[i] = 0;  
  fHistNCMatchvsPt[i] = 0;
  fHistHadCorvsPt[i] = 0;
  fHistNEFvsPtBias[i] = 0;
  fHistNTMatchvsPtBias[i] = 0;  
  fHistNCMatchvsPtBias[i] = 0;
  fHistHadCorvsPtBias[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,"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]);
  }
  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;
    }
    //  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 (centbin == -1)
    centbin = 0;
  if (centbin>2)
    return kTRUE;
  if (!fTracks)
    return kTRUE;

  const Int_t nCluster = fCaloClusters->GetEntriesFast();
  const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
  const Int_t nTrack = fTracks->GetEntriesFast();

  TString fRhoScaledName = fRhoName;
  fRho = GetRhoFromEvent(fRhoScaledName);
  fRhoVal = fRho->GetVal();
  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;
       bool ischecked = false;
       for (Int_t icid = 0;i<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);
     //        fHistNCMatchvsPtvsCent(0),
     fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
     if (jet->MaxTrackPt()<5.0)
       continue;
     fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
     fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
     fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
     if (centbin == 0)
       fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,nTrack);
  }
  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"
	30
	31	ClassImp(AliAnalysisTaskEmcalJetHadCorQA)
	32
	33	//________________________________________________________________________
	34	AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA() :
	35	AliAnalysisTaskEmcalJet("spectra",kFALSE),
	36	fCalo2Name(),
	37	fCaloClusters2(),
	38	fHistRhovsCent(0),
	39	fHistNjetvsCent(0)
	40	{
	41	for (int i = 0;i<3;i++){
	42	fHistNEFvsPt[i] = 0;
	43	fHistNTMatchvsPt[i] = 0;
	44	fHistNCMatchvsPt[i] = 0;
	45	fHistHadCorvsPt[i] = 0;
	46	fHistNEFvsPtBias[i] = 0;
	47	fHistNTMatchvsPtBias[i] = 0;
	48	fHistNCMatchvsPtBias[i] = 0;
	49	fHistHadCorvsPtBias[i] = 0;
	50	}
	51	// Default constructor.
	52
	53	SetMakeGeneralHistograms(kTRUE);
	54	}
	55
	56	//________________________________________________________________________
	57	AliAnalysisTaskEmcalJetHadCorQA::AliAnalysisTaskEmcalJetHadCorQA(const char *name) :
	58	AliAnalysisTaskEmcalJet(name,kTRUE),
	59	fCalo2Name(),
	60	fCaloClusters2(),
	61	fHistRhovsCent(0),
	62	fHistNjetvsCent(0)
	63	{
	64	for (int i = 0;i<3;i++){
65	fHistNEFvsPt[i] = 0;
66	fHistNTMatchvsPt[i] = 0;
67	fHistNCMatchvsPt[i] = 0;
68	fHistHadCorvsPt[i] = 0;
69	fHistNEFvsPtBias[i] = 0;
70	fHistNTMatchvsPtBias[i] = 0;
71	fHistNCMatchvsPtBias[i] = 0;
72	fHistHadCorvsPtBias[i] = 0;
73	}
74	SetMakeGeneralHistograms(kTRUE);
75	}
76
77	//________________________________________________________________________
78	void AliAnalysisTaskEmcalJetHadCorQA::UserCreateOutputObjects()
79	{
80	if (! fCreateHisto)
81	return;
82	AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
83	fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 500, 0, 500);
84	fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
85
86	for (int i = 0;i<3;i++){
87	char name[200];
88	TString title;
89	sprintf(name,"NEFvsPt%i",i);
90	fHistNEFvsPt[i] = new TH2F(name, name, 100,0,1,500,-250,250);
91	fOutput->Add(fHistNEFvsPt[i]);
92	sprintf(name,"NTMatchvsPt%i",i);
93	fHistNTMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
94	fOutput->Add(fHistNTMatchvsPt[i]);
95	sprintf(name,"NCMatchvsPt%i",i);
96	fHistNCMatchvsPt[i] = new TH2F(name, name, 100,0,100,500,-250,250);
97	fOutput->Add(fHistNCMatchvsPt[i]);
98	sprintf(name,"HadCorvsPt%i",i);
99	fHistHadCorvsPt[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
100	fOutput->Add(fHistHadCorvsPt[i]);
101	sprintf(name,"NEFvsPtBias%i",i);
102	fHistNEFvsPtBias[i] = new TH2F(name, name, 100,0,1,500,-250,250);
103	fOutput->Add(fHistNEFvsPtBias[i]);
104	sprintf(name,"NTMatchvsPtBias%i",i);
105	fHistNTMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
106	fOutput->Add(fHistNTMatchvsPtBias[i]);
107	sprintf(name,"NCMatchvsPtBias%i",i);
108	fHistNCMatchvsPtBias[i] = new TH2F(name, name, 100,0,100,500,-250,250);
109	fOutput->Add(fHistNCMatchvsPtBias[i]);
110	sprintf(name,"HadCorvsPtBias%i",i);
111	fHistHadCorvsPtBias[i] = new TH2F(name, name, 1000,0,500,500,-250,250);
112	fOutput->Add(fHistHadCorvsPtBias[i]);
113	}
114	fHistNTMatchvsPtvsNtack0 = new TH3F("NTMmatchvsPtvsNtrack0", "NTMatchsvsPtvsNtrack0", 100,0,100,500,-250,250,250,0,2500);
115
116	// for (Int_t i = 0;i<6;++i){
117	// name = TString(Form("JetPtvsTrackPt_%i",i));
118	// title = TString(Form("Jet pT vs Leading Track pT cent bin %i",i));
119	// fHistJetPtvsTrackPt[i] = new TH2F(name,title,1000,-500,500,100,0,100);
120	// fOutput->Add(fHistJetPtvsTrackPt[i]);
121
122	// }
123	fOutput->Add(fHistRhovsCent);
124	fOutput->Add(fHistNjetvsCent);
125	fOutput->Add(fHistNTMatchvsPtvsNtack0);
126	PostData(1, fOutput);
127	}
128
129	//________________________________________________________________________
130
131	Int_t AliAnalysisTaskEmcalJetHadCorQA::GetCentBin(Double_t cent) const
132	{
133	// Get centrality bin.
134
135	Int_t centbin = -1;
136	if (cent>=0 && cent<10)
137	centbin = 0;
138	else if (cent>=10 && cent<30)
139	centbin = 1;
140	else if (cent>=30 && cent<50)
141	centbin = 2;
142	else if (cent>50)
143	centbin =3;
144	return centbin;
145	}
146
147	//________________________________________________________________________
148
149	Float_t AliAnalysisTaskEmcalJetHadCorQA:: RelativePhi(Double_t mphi,Double_t vphi) const
150	{
151	if (vphi < -1TMath::Pi()) vphi += (2TMath::Pi());
152	else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
153	if (mphi < -1TMath::Pi()) mphi += (2TMath::Pi());
154	else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
155	double dphi = mphi-vphi;
156	if (dphi < -1TMath::Pi()) dphi += (2TMath::Pi());
157	else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
158
159	return dphi;//dphi in [-Pi, Pi]
160	}
161
162	//________________________________________________________________________
163	void AliAnalysisTaskEmcalJetHadCorQA::ExecOnce(){
164	AliAnalysisTaskEmcalJet::ExecOnce();
165
166	// AliAnalysisManger *am = AliAnalysisManager::GetAnalysisManger();
167	// if (fCaloName == "CaloClusters")
168	// am->LoadBranch("CaloClusters");
169
170	if (!fCalo2Name.IsNull() && !fCaloClusters2){
171	fCaloClusters2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
172	if (!fCaloClusters2){
173	AliError(Form("%s: Could not retrieve calo clusters %s!",GetName(),fCalo2Name.Data()));
174	fInitialized = kFALSE;
175	return;
176	}
177	// fCaloClusters2(),
178	// else if (!fJets->GetClass()->GetBaseClass("AliVCluster")){
179	// AliError(Form("%s: Collection %s does not contain AliEmcalParticle objects!",GetName(),fCalo2Name.Data()));
180	// fCaloClusters2 = 0;
181	// fInitialized = kFALSE;
182	// return;
183	// }
184	}
185	}
186
187
188	//________________________________________________________________________
189	Bool_t AliAnalysisTaskEmcalJetHadCorQA::Run()
190	{
191	Int_t centbin = GetCentBin(fCent);
192	//for pp analyses we will just use the first centrality bin
193	if (centbin == -1)
194	centbin = 0;
195	if (centbin>2)
196	return kTRUE;
197	if (!fTracks)
198	return kTRUE;
199
200	const Int_t nCluster = fCaloClusters->GetEntriesFast();
201	const Int_t nCluster2 = fCaloClusters2->GetEntriesFast();
202	const Int_t nTrack = fTracks->GetEntriesFast();
203
204	TString fRhoScaledName = fRhoName;
205	fRho = GetRhoFromEvent(fRhoScaledName);
206	fRhoVal = fRho->GetVal();
207	fHistRhovsCent->Fill(fCent,fRhoVal);
208	const Int_t Njets = fJets->GetEntriesFast();
209
210	Int_t NjetAcc = 0;
211
212
213	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
214	Int_t TrackMatch = 0;
215	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(iJets));
216	if (!jet)
217	continue;
218	if (jet->Area()==0)
219	continue;
220	if (jet->Pt()<0.1)
221	continue;
222	if (jet->MaxTrackPt()>100)
223	continue;
224	if (! AcceptJet(jet))
225	continue;
226	NjetAcc++;
227	vector<Int_t> cluster_id; //we need to keep track of the jet clusters that we find
228	vector<Int_t> cluster_id2;
229	Double_t Esub = 0; //total E subtracted from the jet
230	Double_t jetPt = -500;
231	jetPt = jet->Pt()-jet->Area()*fRhoVal;
232	for (int i = 0;i<nTrack;i++){
233	AliEmcalParticle emctrack = static_cast<AliEmcalParticle>(fTracks->At(i));
234	if (!emctrack)
235	continue;
236	if (!emctrack->IsEMCAL())
237	continue;
238	AliVTrack track = (AliVTrack)emctrack->GetTrack();
239	if (! track)
240	continue;
241	if (! AcceptTrack(track))
242	continue;
243	if (! IsJetTrack(jet,i,false))
244	continue;
245	Int_t iClus = track->GetEMCALcluster();
246	if (iClus<0)
247	continue;
248	bool ischecked = false;
249	for (Int_t icid = 0;i<cluster_id.size();icid++)
250	if (cluster_id[icid] == iClus)
251	ischecked = true; // we've already looked at this uncorrected cluster
252	if (ischecked)
253	continue; //no need to go further
254	AliEmcalParticle emcluster = static_cast<AliEmcalParticle>(fCaloClusters->At(iClus));
255	AliVCluster *cluster = emcluster->GetCluster();
256	if (! cluster)
257	continue;
258	if (! AcceptCluster(cluster))
259	continue;
260	Double_t etadiff = 999;
261	Double_t phidiff = 999;
262	AliPicoTrack::GetEtaPhiDiff(track,cluster,phidiff,etadiff);
263	if (! (TMath::Abs(phidiff)< 0.025&&TMath::Abs(etadiff)<0.015))
264	continue;
265	TrackMatch++; //this cluster has been matched, let's add it to the list
266	cluster_id.push_back(iClus);
267	Int_t ismatch = -1;
268	//now we need to find its matched corrected cluster if any
269	for (Int_t ic = 0;ic < nCluster2;ic++){
270	//we don't need to check the list of corrected clusters because 1 to 1 between corrected and uncorrected
271	AliVCluster emcluster2 = static_cast<AliVCluster>(fCaloClusters2->At(ic));
272	if (! emcluster2)
273	continue;
274	if (! AcceptCluster(emcluster2))
275	continue;
276	if (! IsJetCluster(jet,ic,false))
277	continue;
278	TLorentzVector nPart;
279	cluster->GetMomentum(nPart,const_cast<Double_t*>(fVertex));
280	TLorentzVector nPart2;
281	emcluster2->GetMomentum(nPart2,const_cast<Double_t*>(fVertex));
282	float R = pow(pow(nPart.Eta()-nPart2.Eta(),2)+pow(nPart.Phi()-nPart2.Phi(),2),0.5);
283	if (R < 0.001){//this cluster stayed in the jet!
284	cluster_id2.push_back(ic);
285	ismatch = ic;
286	}
287	} // end of cluster2 loop
288	//we only get here if the track was matched to a cluster that hadn't been looked at before
289	if (ismatch < 0) //this cluster was entirely deleted
290	Esub+=cluster->E();
291	else{ //get the corrected cluster
292	AliVCluster emcluster2temp = static_cast<AliVCluster>(fCaloClusters2->At(ismatch));
293	Esub+=(cluster->E() - emcluster2temp->E());
294	}
295
296	} // end of track loop
297	fHistNEFvsPt[centbin]->Fill(jet->NEF(),jetPt);
298	fHistNTMatchvsPt[centbin]->Fill(TrackMatch,jetPt);
299	// fHistNCMatchvsPtvsCent(0),
300	fHistHadCorvsPt[centbin]->Fill(Esub,jetPt);
301	if (jet->MaxTrackPt()<5.0)
302	continue;
303	fHistNEFvsPtBias[centbin]->Fill(jet->NEF(),jetPt);
304	fHistNTMatchvsPtBias[centbin]->Fill(TrackMatch,jetPt);
305	fHistHadCorvsPtBias[centbin]->Fill(Esub,jetPt);
306	if (centbin == 0)
307	fHistNTMatchvsPtvsNtack0->Fill(TrackMatch,jetPt,nTrack);
308	}
309	fHistNjetvsCent->Fill(fCent,NjetAcc);
310	return kTRUE;
311	}
312
313
314
315
316