[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliHadCorrTask.cxx

// $Id$
//
// Hadronic correction task.
//
//

#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TParticle.h>
#include "AliAODCaloCluster.h"
#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
#include "AliCentrality.h"
#include "AliESDCaloCluster.h"
#include "AliESDtrack.h"
#include "AliFJWrapper.h"
#include "AliPicoTrack.h"
#include "AliVEventHandler.h"
#include "AliHadCorrTask.h"

ClassImp(AliHadCorrTask)

//________________________________________________________________________
AliHadCorrTask::AliHadCorrTask() : 
  AliAnalysisTaskSE("AliHadCorrTask"),
  fTracksName(),
  fCaloName(),
  fOutCaloName(),
  fHadCorr(0),
  fMinPt(0.15),
  fOutClusters(0),
  fOutputList(0),
  fHistNclusvsCent(0),
  fHistNclusMatchvsCent(0),
  fHistEbefore(0),
  fHistEafter(0),
  fHistEoPCent(0),
  fHistNMatchCent(0)
{
  // Default constructor.

  for(Int_t i=0; i<4; ++i) {
    fHistMatchEvsP[i]   = 0;
    fHistMatchdRvsEP[i] = 0;
    for(Int_t j=0; j<5; ++j) 
      fHistMatchEtaPhi[i][j] = 0;
  }
}

//________________________________________________________________________
AliHadCorrTask::AliHadCorrTask(const char *name) : 
  AliAnalysisTaskSE(name),
  fTracksName("Tracks"),
  fCaloName("CaloClusters"),
  fOutCaloName("CaloClustersCorr"),
  fHadCorr(0),
  fMinPt(0.15),
  fOutClusters(0),
  fOutputList(0),
  fHistNclusvsCent(0),
  fHistNclusMatchvsCent(0),
  fHistEbefore(0),
  fHistEafter(0),
  fHistEoPCent(0),
  fHistNMatchCent(0)
{
  // Standard constructor.

  for(Int_t i=0; i<4; ++i) {
    fHistMatchEvsP[i]   = 0;
    fHistMatchdRvsEP[i] = 0;
    for(Int_t j=0; j<5; ++j) 
      fHistMatchEtaPhi[i][j] = 0;
  }

  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";

  DefineInput(0,TChain::Class());
  DefineOutput(1,TList::Class());
}

//________________________________________________________________________
AliHadCorrTask::~AliHadCorrTask()
{
  // Destructor
}

//________________________________________________________________________
Int_t AliHadCorrTask::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 && cent<=100) 
    centbin=3;

  return centbin;
}

//________________________________________________________________________
Int_t AliHadCorrTask::GetMomBin(Double_t p) const
{
  // Get momenum bin.

  Int_t pbin=-1;
  if (p>=0 && p<0.5) 
    pbin=0;
  else if (p>=0.5 && p<2.) 
    pbin=1;
  else if (p>=2. && p<3.) 
    pbin=2;
  else if (p>=3. && p<5.) 
    pbin=3;
  else if (p>=5.) 
    pbin=4;

  return pbin;
}

//________________________________________________________________________
void AliHadCorrTask::UserCreateOutputObjects()
{
  // Create my user objects.

  AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
  if (!handler) {
    AliError("Input handler not available!");
    return;
  }
 
  if (handler->InheritsFrom("AliESDInputHandler")) {
    fOutClusters = new TClonesArray("AliESDCaloCluster");
  } else if (handler->InheritsFrom("AliAODInputHandler")) {
    fOutClusters = new TClonesArray("AliAODCaloCluster");
  } else {
    AliError("Input handler not recognized!");
    return;
  }
  fOutClusters->SetName(fOutCaloName);
 
  OpenFile(1);
  fOutputList = new TList();
  fOutputList->SetOwner();

  for(Int_t icent=0; icent<4; ++icent) {
    for(Int_t ipt=0; ipt<5; ++ipt){
      TString name(Form("fHistMatchEtaPhi_%i_%i",icent,ipt));
      fHistMatchEtaPhi[icent][ipt] = new TH2F(name, name, 400, -0.2, 0.2, 1600, -0.8, 0.8);
      fOutputList->Add(fHistMatchEtaPhi[icent][ipt]);
    }

    TString name(Form("fHistMatchEvsP_%i",icent));
    fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
    fOutputList->Add(fHistMatchEvsP[icent]);

    name = Form("fHistMatchdRvsEP_%i",icent);
    fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
    fOutputList->Add(fHistMatchdRvsEP[icent]);
  }

  fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
  fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
  fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
  fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
  fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 101, -0.5, 100.5);
  fOutputList->Add(fHistNclusMatchvsCent);
  fOutputList->Add(fHistNclusvsCent);
  fOutputList->Add(fHistEbefore);
  fOutputList->Add(fHistEafter);
  fOutputList->Add(fHistEoPCent);
  fOutputList->Add(fHistNMatchCent);

  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliHadCorrTask::UserExec(Option_t *) 
{
  // Execute per event.

  // post output in event if not yet present
  if (!(InputEvent()->FindListObject(fOutCaloName)))
    InputEvent()->AddObject(fOutClusters);
  
  // delete output
  fOutClusters->Delete();

  // esd or aod mode
  Bool_t esdMode = kTRUE;
  if (dynamic_cast<AliAODEvent*>(InputEvent()))
      esdMode = kFALSE;

  // optimization in case autobranch loading is off
  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  if (fCaloName == "CaloClusters")
    am->LoadBranch("CaloClusters");
  if (fTracksName == "Tracks")
    am->LoadBranch("Tracks");

  // get centrality 
  Double_t cent = -1; 
  AliCentrality *centrality = InputEvent()->GetCentrality() ;
  if (centrality)
    cent = centrality->GetCentralityPercentile("V0M");
  else
    cent=99; // probably pp data
  if (cent<0) {
    AliError(Form("Centrality negative: %f", cent));
    return;
  }
  Int_t centbin = GetCentBin(cent);

  // get input collections
  TClonesArray *tracks = 0;
  TClonesArray *clus   = 0;
  TList *l = InputEvent()->GetList();
  tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
  if (!tracks) {
    AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
    return;
  }
  clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
  if (!clus) {
    AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
    return;
  }

  // get primary vertex
  Double_t vertex[3] = {0, 0, 0};
  InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);

  // loop over clusters and tracks
  const Int_t Nclus = clus->GetEntries();
  const Int_t Ntrks = tracks->GetEntries();
  for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
    AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
    if (!c)
      continue;
    if (!c->IsEMCAL())
      continue;

    // make primary particle out of cluster
    TLorentzVector nPart;
    c->GetMomentum(nPart, vertex);

    Double_t etclus = nPart.Pt();
    if (etclus<fMinPt) 
      continue;

    Double_t energyclus = nPart.P();

    // reset cluster/track matching
    c->SetEmcCpvDistance(-1);
    c->SetTrackDistance(999,999);

    // run cluster-track matching
    Int_t    imin       = -1;
    Double_t dEtaMin    = 1e9;
    Double_t dPhiMin    = 1e9;
    Double_t dRmin      = 1e9;
    Double_t totalTrkP  = 0.0; // count total track momentum
    Int_t    Nmatches   = 0;   // count total number of matches
    for (Int_t t = 0; t<Ntrks; ++t) {
      AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(t));
      if (!track)
        continue;
      if (!track->IsEMCAL())
        continue;
      if (track->Pt()<fMinPt)
        continue;

      Double_t etadiff=999;
      Double_t phidiff=999;
      AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
      Double_t dR = TMath::Sqrt(etadiff*etadiff+phidiff*phidiff);
      if(dR<dRmin){
        dEtaMin = etadiff;
        dPhiMin = phidiff;
        dRmin = dR;
        imin = t;
      }
      if (fHadCorr>1) {
        Double_t mom = track->P();
        Int_t mombin = GetMomBin(mom);
        if (mombin>-1) {
          fHistMatchEtaPhi[centbin][mombin]->Fill(etadiff,phidiff);
          fHistMatchdRvsEP[centbin]->Fill(dR,energyclus/mom);
        }
      }
      if (TMath::Abs(phidiff)<0.05 && TMath::Abs(etadiff)<0.025) { // pp cuts!!!
        ++Nmatches;
        totalTrkP += track->P();
      }
    }

    // store closest track
    c->SetEmcCpvDistance(imin);
    c->SetTrackDistance(dPhiMin, dEtaMin);

    fHistNclusvsCent->Fill(cent);
    fHistEbefore->Fill(cent,energyclus);
    fHistNMatchCent->Fill(cent,Nmatches);
    if(Nmatches>0) 
      fHistNclusMatchvsCent->Fill(cent);

    // apply correction / subtraction
    if (fHadCorr>0) {
      // to subtract only the closest track set fHadCor to a %
      // to subtract all tracks within the cut set fHadCor to %+1
      if (fHadCorr>1) {
        if (totalTrkP>0) {
          Double_t EoP = energyclus/totalTrkP;
          fHistEoPCent->Fill(cent,EoP);
          fHistMatchEvsP[centbin]->Fill(energyclus,EoP);
        }
        energyclus -= (fHadCorr-1)*totalTrkP;
      } else if (imin>=0) {
        AliVTrack *t = dynamic_cast<AliVTrack*>(tracks->At(imin));
        if (t) {
          Double_t mom = t->P();
          Int_t mombin = GetMomBin(mom);
          fHistMatchEtaPhi[centbin][mombin]->Fill(dEtaMin,dPhiMin);
          if (mom>0){
            fHistMatchEvsP[centbin]->Fill(energyclus,energyclus/mom);
            fHistEoPCent->Fill(cent,energyclus/mom);
            fHistMatchdRvsEP[centbin]->Fill(dRmin,energyclus/mom);
          }
          if (TMath::Abs(dPhiMin)<0.05 && TMath::Abs(dEtaMin)<0.025) { // pp cuts!!!
            energyclus -= fHadCorr*t->P();
          }
        }
      }
    }
    if (energyclus<0) 
      energyclus = 0;
    fHistEafter->Fill(cent,energyclus);

    if (energyclus>0) { // create corrected cluster
      AliVCluster *oc;
      if (esdMode) {
        oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*>(c)));
      } else {
        oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*>(c)));
      }
      oc->SetE(energyclus);
      ++clusCount;
    }
  }
}

//________________________________________________________________________
void AliHadCorrTask::Terminate(Option_t *) 
{
  // Nothing to be done.
}
Commit	Line	Data
d6f334b5	1	// $Id$
	2	//
	3	// Hadronic correction task.
	4	//
	5	//
	6
	7	#include <TChain.h>
0b777a09	8	#include <TClonesArray.h>
35789a2d	9	#include <TH1F.h>
35789a2d	10	#include <TH2F.h>
d6f334b5	11	#include <TList.h>
35789a2d	12	#include <TLorentzVector.h>
d6f334b5	13	#include <TParticle.h>
	14	#include "AliAODCaloCluster.h"
	15	#include "AliAODEvent.h"
0b777a09	16	#include "AliAnalysisManager.h"
d6f334b5	17	#include "AliCentrality.h"
d6f334b5	18	#include "AliESDCaloCluster.h"
0b777a09	19	#include "AliESDtrack.h"
0b777a09	20	#include "AliFJWrapper.h"
0b777a09	21	#include "AliPicoTrack.h"
35789a2d	22	#include "AliVEventHandler.h"
35789a2d	23	#include "AliHadCorrTask.h"
0b777a09	24
	25	ClassImp(AliHadCorrTask)
	26
	27	//________________________________________________________________________
d6f334b5	28	AliHadCorrTask::AliHadCorrTask() :
0b777a09	29	AliAnalysisTaskSE("AliHadCorrTask"),
d6f334b5	30	fTracksName(),
	31	fCaloName(),
	32	fOutCaloName(),
f22cb876	33	fHadCorr(0),
0b777a09	34	fMinPt(0.15),
f22cb876	35	fOutClusters(0),
0b777a09	36	fOutputList(0),
0b777a09	37	fHistNclusvsCent(0),
0b777a09	38	fHistNclusMatchvsCent(0),
f22cb876	39	fHistEbefore(0),
4711c521	40	fHistEafter(0),
	41	fHistEoPCent(0),
	42	fHistNMatchCent(0)
0b777a09	43	{
d6f334b5	44	// Default constructor.
0b777a09	45
d6f334b5	46	for(Int_t i=0; i<4; ++i) {
	47	fHistMatchEvsP[i] = 0;
	48	fHistMatchdRvsEP[i] = 0;
	49	for(Int_t j=0; j<5; ++j)
	50	fHistMatchEtaPhi[i][j] = 0;
0b777a09	51	}
0b777a09	52	}
0b777a09	53
d6f334b5	54	//________________________________________________________________________
	55	AliHadCorrTask::AliHadCorrTask(const char *name) :
	56	AliAnalysisTaskSE(name),
0b777a09	57	fTracksName("Tracks"),
0b777a09	58	fCaloName("CaloClusters"),
d6f334b5	59	fOutCaloName("CaloClustersCorr"),
f22cb876	60	fHadCorr(0),
f22cb876	61	fMinPt(0.15),
f22cb876	62	fOutClusters(0),
	63	fOutputList(0),
	64	fHistNclusvsCent(0),
	65	fHistNclusMatchvsCent(0),
	66	fHistEbefore(0),
4711c521	67	fHistEafter(0),
	68	fHistEoPCent(0),
	69	fHistNMatchCent(0)
0b777a09	70	{
	71	// Standard constructor.
	72
d6f334b5	73	for(Int_t i=0; i<4; ++i) {
	74	fHistMatchEvsP[i] = 0;
	75	fHistMatchdRvsEP[i] = 0;
	76	for(Int_t j=0; j<5; ++j)
	77	fHistMatchEtaPhi[i][j] = 0;
	78	}
	79
0b777a09	80	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
0b777a09	81
d6f334b5	82	DefineInput(0,TChain::Class());
	83	DefineOutput(1,TList::Class());
	84	}
0b777a09	85
	86	//________________________________________________________________________
	87	AliHadCorrTask::~AliHadCorrTask()
	88	{
	89	// Destructor
	90	}
	91
d6f334b5	92	//________________________________________________________________________
	93	Int_t AliHadCorrTask::GetCentBin(Double_t cent) const
	94	{
	95	// Get centrality bin.
	96
	97	Int_t centbin = -1;
	98	if (cent>=0 && cent<10)
	99	centbin=0;
	100	else if (cent>=10 && cent<30)
	101	centbin=1;
	102	else if (cent>=30 && cent<50)
	103	centbin=2;
	104	else if (cent>=50 && cent<=100)
	105	centbin=3;
	106
	107	return centbin;
	108	}
	109
	110	//________________________________________________________________________
	111	Int_t AliHadCorrTask::GetMomBin(Double_t p) const
	112	{
	113	// Get momenum bin.
	114
	115	Int_t pbin=-1;
	116	if (p>=0 && p<0.5)
	117	pbin=0;
	118	else if (p>=0.5 && p<2.)
	119	pbin=1;
	120	else if (p>=2. && p<3.)
	121	pbin=2;
	122	else if (p>=3. && p<5.)
	123	pbin=3;
	124	else if (p>=5.)
	125	pbin=4;
	126
	127	return pbin;
	128	}
	129
0b777a09	130	//________________________________________________________________________
	131	void AliHadCorrTask::UserCreateOutputObjects()
	132	{
d6f334b5	133	// Create my user objects.
0b777a09	134
35789a2d	135	AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
35789a2d	136	if (!handler) {
	137	AliError("Input handler not available!");
	138	return;
	139	}
	140
	141	if (handler->InheritsFrom("AliESDInputHandler")) {
	142	fOutClusters = new TClonesArray("AliESDCaloCluster");
d6f334b5	143	} else if (handler->InheritsFrom("AliAODInputHandler")) {
35789a2d	144	fOutClusters = new TClonesArray("AliAODCaloCluster");
d6f334b5	145	} else {
35789a2d	146	AliError("Input handler not recognized!");
	147	return;
	148	}
0b777a09	149	fOutClusters->SetName(fOutCaloName);
	150
	151	OpenFile(1);
	152	fOutputList = new TList();
	153	fOutputList->SetOwner();
	154
d6f334b5	155	for(Int_t icent=0; icent<4; ++icent) {
	156	for(Int_t ipt=0; ipt<5; ++ipt){
	157	TString name(Form("fHistMatchEtaPhi_%i_%i",icent,ipt));
	158	fHistMatchEtaPhi[icent][ipt] = new TH2F(name, name, 400, -0.2, 0.2, 1600, -0.8, 0.8);
4711c521	159	fOutputList->Add(fHistMatchEtaPhi[icent][ipt]);
	160	}
	161
d6f334b5	162	TString name(Form("fHistMatchEvsP_%i",icent));
d6f334b5	163	fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
4711c521	164	fOutputList->Add(fHistMatchEvsP[icent]);
0b777a09	165
d6f334b5	166	name = Form("fHistMatchdRvsEP_%i",icent);
	167	fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
	168	fOutputList->Add(fHistMatchdRvsEP[icent]);
	169	}
0b777a09	170
d6f334b5	171	fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent", 100, 0, 100);
	172	fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
	173	fHistEbefore = new TH1F("Ebefore", "Ebefore", 100, 0, 100);
	174	fHistEafter = new TH1F("Eafter", "Eafter", 100, 0, 100);
	175	fHistEoPCent = new TH2F("EoPCent", "EoPCent", 100, 0, 100, 1000, 0, 10);
	176	fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent", 100, 0, 100, 101, -0.5, 100.5);
0b777a09	177	fOutputList->Add(fHistNclusMatchvsCent);
	178	fOutputList->Add(fHistNclusvsCent);
	179	fOutputList->Add(fHistEbefore);
	180	fOutputList->Add(fHistEafter);
4711c521	181	fOutputList->Add(fHistEoPCent);
4711c521	182	fOutputList->Add(fHistNMatchCent);
0b777a09	183
	184	PostData(1, fOutputList);
	185	}
	186
	187	//________________________________________________________________________
	188	void AliHadCorrTask::UserExec(Option_t *)
	189	{
d6f334b5	190	// Execute per event.
0b777a09	191
d6f334b5	192	// post output in event if not yet present
0b777a09	193	if (!(InputEvent()->FindListObject(fOutCaloName)))
0b777a09	194	InputEvent()->AddObject(fOutClusters);
d6f334b5	195
	196	// delete output
	197	fOutClusters->Delete();
0b777a09	198
d6f334b5	199	// esd or aod mode
	200	Bool_t esdMode = kTRUE;
	201	if (dynamic_cast<AliAODEvent*>(InputEvent()))
	202	esdMode = kFALSE;
0b777a09	203
d6f334b5	204	// optimization in case autobranch loading is off
0b777a09	205	AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
d6f334b5	206	if (fCaloName == "CaloClusters")
	207	am->LoadBranch("CaloClusters");
	208	if (fTracksName == "Tracks")
	209	am->LoadBranch("Tracks");
0b777a09	210
d6f334b5	211	// get centrality
c7b78119	212	Double_t cent = -1;
d6f334b5	213	AliCentrality *centrality = InputEvent()->GetCentrality() ;
	214	if (centrality)
	215	cent = centrality->GetCentralityPercentile("V0M");
0b777a09	216	else
d6f334b5	217	cent=99; // probably pp data
	218	if (cent<0) {
	219	AliError(Form("Centrality negative: %f", cent));
	220	return;
	221	}
	222	Int_t centbin = GetCentBin(cent);
0b777a09	223
d6f334b5	224	// get input collections
	225	TClonesArray *tracks = 0;
	226	TClonesArray *clus = 0;
	227	TList *l = InputEvent()->GetList();
0b777a09	228	tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
	229	if (!tracks) {
	230	AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
	231	return;
	232	}
0b777a09	233	clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
	234	if (!clus) {
	235	AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
	236	return;
	237	}
	238
d6f334b5	239	// get primary vertex
	240	Double_t vertex[3] = {0, 0, 0};
	241	InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
	242
	243	// loop over clusters and tracks
	244	const Int_t Nclus = clus->GetEntries();
	245	const Int_t Ntrks = tracks->GetEntries();
	246	for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
	247	AliVCluster c = dynamic_cast<AliVCluster>(clus->At(iClus));
	248	if (!c)
	249	continue;
	250	if (!c->IsEMCAL())
	251	continue;
	252
	253	// make primary particle out of cluster
	254	TLorentzVector nPart;
	255	c->GetMomentum(nPart, vertex);
	256
	257	Double_t etclus = nPart.Pt();
	258	if (etclus<fMinPt)
	259	continue;
	260
	261	Double_t energyclus = nPart.P();
	262
	263	// reset cluster/track matching
	264	c->SetEmcCpvDistance(-1);
	265	c->SetTrackDistance(999,999);
	266
	267	// run cluster-track matching
	268	Int_t imin = -1;
	269	Double_t dEtaMin = 1e9;
	270	Double_t dPhiMin = 1e9;
cce6a19c	271	Double_t dRmin = 1e9;
d6f334b5	272	Double_t totalTrkP = 0.0; // count total track momentum
	273	Int_t Nmatches = 0; // count total number of matches
	274	for (Int_t t = 0; t<Ntrks; ++t) {
	275	AliVTrack track = dynamic_cast<AliVTrack>(tracks->At(t));
	276	if (!track)
0b777a09	277	continue;
d6f334b5	278	if (!track->IsEMCAL())
0b777a09	279	continue;
d6f334b5	280	if (track->Pt()<fMinPt)
	281	continue;
	282
	283	Double_t etadiff=999;
	284	Double_t phidiff=999;
	285	AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
	286	Double_t dR = TMath::Sqrt(etadiffetadiff+phidiffphidiff);
	287	if(dR<dRmin){
	288	dEtaMin = etadiff;
	289	dPhiMin = phidiff;
	290	dRmin = dR;
	291	imin = t;
	292	}
	293	if (fHadCorr>1) {
	294	Double_t mom = track->P();
	295	Int_t mombin = GetMomBin(mom);
	296	if (mombin>-1) {
	297	fHistMatchEtaPhi[centbin][mombin]->Fill(etadiff,phidiff);
	298	fHistMatchdRvsEP[centbin]->Fill(dR,energyclus/mom);
	299	}
0b777a09	300	}
d6f334b5	301	if (TMath::Abs(phidiff)<0.05 && TMath::Abs(etadiff)<0.025) { // pp cuts!!!
	302	++Nmatches;
	303	totalTrkP += track->P();
	304	}
	305	}
	306
	307	// store closest track
	308	c->SetEmcCpvDistance(imin);
	309	c->SetTrackDistance(dPhiMin, dEtaMin);
	310
	311	fHistNclusvsCent->Fill(cent);
	312	fHistEbefore->Fill(cent,energyclus);
	313	fHistNMatchCent->Fill(cent,Nmatches);
	314	if(Nmatches>0)
	315	fHistNclusMatchvsCent->Fill(cent);
	316
	317	// apply correction / subtraction
	318	if (fHadCorr>0) {
	319	// to subtract only the closest track set fHadCor to a %
	320	// to subtract all tracks within the cut set fHadCor to %+1
	321	if (fHadCorr>1) {
	322	if (totalTrkP>0) {
	323	Double_t EoP = energyclus/totalTrkP;
	324	fHistEoPCent->Fill(cent,EoP);
	325	fHistMatchEvsP[centbin]->Fill(energyclus,EoP);
	326	}
	327	energyclus -= (fHadCorr-1)*totalTrkP;
	328	} else if (imin>=0) {
	329	AliVTrack t = dynamic_cast<AliVTrack>(tracks->At(imin));
	330	if (t) {
	331	Double_t mom = t->P();
	332	Int_t mombin = GetMomBin(mom);
	333	fHistMatchEtaPhi[centbin][mombin]->Fill(dEtaMin,dPhiMin);
	334	if (mom>0){
	335	fHistMatchEvsP[centbin]->Fill(energyclus,energyclus/mom);
	336	fHistEoPCent->Fill(cent,energyclus/mom);
	337	fHistMatchdRvsEP[centbin]->Fill(dRmin,energyclus/mom);
	338	}
	339	if (TMath::Abs(dPhiMin)<0.05 && TMath::Abs(dEtaMin)<0.025) { // pp cuts!!!
	340	energyclus -= fHadCorr*t->P();
	341	}
	342	}
0b777a09	343	}
0b777a09	344	}
d6f334b5	345	if (energyclus<0)
	346	energyclus = 0;
	347	fHistEafter->Fill(cent,energyclus);
	348
	349	if (energyclus>0) { // create corrected cluster
	350	AliVCluster *oc;
	351	if (esdMode) {
	352	oc = new ((fOutClusters)[clusCount]) AliESDCaloCluster((dynamic_cast<AliESDCaloCluster*>(c)));
	353	} else {
	354	oc = new ((fOutClusters)[clusCount]) AliAODCaloCluster((dynamic_cast<AliAODCaloCluster*>(c)));
	355	}
	356	oc->SetE(energyclus);
	357	++clusCount;
	358	}
0b777a09	359	}
0b777a09	360	}
4711c521	361
0b777a09	362	//________________________________________________________________________
	363	void AliHadCorrTask::Terminate(Option_t *)
	364	{
d6f334b5	365	// Nothing to be done.
0b777a09	366	}
0b777a09	367