//
// Hadronic correction task.
//
-// Author: R.Reed, C.Loizides
+// Author: R.Reed, C.Loizides, S.Aiola
#include <TChain.h>
#include <TClonesArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TList.h>
-#include <TLorentzVector.h>
-#include "AliAODCaloCluster.h"
-#include "AliAODEvent.h"
#include "AliAnalysisManager.h"
-#include "AliCentrality.h"
+#include "AliAODEvent.h"
+#include "AliAODCaloCluster.h"
#include "AliESDCaloCluster.h"
-#include "AliESDtrack.h"
+#include "AliVTrack.h"
#include "AliPicoTrack.h"
#include "AliVEventHandler.h"
+#include "AliEmcalParticle.h"
#include "AliHadCorrTask.h"
//________________________________________________________________________
AliHadCorrTask::AliHadCorrTask() :
- AliAnalysisTaskSE("AliHadCorrTask"),
- fTracksName(),
- fCaloName(),
+ AliAnalysisTaskEmcal("AliHadCorrTask", kFALSE),
fOutCaloName(),
fPhiMatch(0.05),
fEtaMatch(0.025),
fDoTrackClus(0),
fHadCorr(0),
- fMinPt(0.15),
- fCreateHisto(kFALSE),
+ fEexclCell(0),
fOutClusters(0),
- fOutputList(0),
fHistNclusvsCent(0),
fHistNclusMatchvsCent(0),
fHistEbefore(0),
fHistEafter(0),
fHistEoPCent(0),
fHistNMatchCent(0),
- fHistNMatchCent_trk(0),
fHistCentrality(0)
{
// Default constructor.
for(Int_t i=0; i<8; i++) {
+ fHistEsubPch[i] = 0;
+ fHistEsubPchRat[i] = 0;
+ for(Int_t j=0; j<4; j++) {
+ fHistNCellsEnergy[i][j] = 0;
+ }
if (i<4) {
- fHistMatchEvsP[i] = 0;
- fHistMatchdRvsEP[i] = 0;
- for(Int_t j=0; j<3; j++) {
- fHistEsubPch[i][j] = 0;
- fHistEsubPchRat[i][j] = 0;
- }
+ fHistMatchEvsP[i] = 0;
+ fHistMatchdRvsEP[i] = 0;
+ fHistNMatchEnergy[i] = 0;
}
for(Int_t j=0; j<9; j++) {
- fHistMatchEtaPhi[i][j] = 0;
+ for(Int_t k=0; k<2; k++) {
+ fHistMatchEtaPhi[i][j][k] = 0;
+ }
}
-
}
}
//________________________________________________________________________
-AliHadCorrTask::AliHadCorrTask(const char *name) :
- AliAnalysisTaskSE(name),
- fTracksName("Tracks"),
- fCaloName("CaloClusters"),
+AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
+ AliAnalysisTaskEmcal(name, histo),
fOutCaloName("CaloClustersCorr"),
fPhiMatch(0.05),
fEtaMatch(0.025),
- fDoTrackClus(0),
+ fDoTrackClus(1),
fHadCorr(0),
- fMinPt(0.15),
- fCreateHisto(kFALSE),
+ fEexclCell(0),
fOutClusters(0),
- fOutputList(0),
fHistNclusvsCent(0),
fHistNclusMatchvsCent(0),
fHistEbefore(0),
fHistEafter(0),
fHistEoPCent(0),
fHistNMatchCent(0),
- fHistNMatchCent_trk(0),
fHistCentrality(0)
-
{
// Standard constructor.
fHistMatchEvsP[i] = 0;
fHistMatchdRvsEP[i] = 0;
for(Int_t j=0; j<3; j++) {
- fHistEsubPch[i][j] = 0;
- fHistEsubPchRat[i][j] = 0;
+ fHistNCellsEnergy[i][j] = 0;
}
}
+ fHistEsubPch[i] = 0;
+ fHistEsubPchRat[i] = 0;
for(Int_t j=0; j<9; j++) {
- fHistMatchEtaPhi[i][j] = 0;
+ for(Int_t k=0; k<2; k++) {
+ fHistMatchEtaPhi[i][j][k] = 0;
+ }
}
}
fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
-
- DefineOutput(1, TList::Class());
}
//________________________________________________________________________
// 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)
+ if (p<0.5)
pbin=0;
else if (p>=0.5 && p<1.0)
pbin=1;
//________________________________________________________________________
Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
{
+ // Get sigma in eta.
- Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.045,0.042};
+ Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
return 2.0*EtaSigma[pbin];
}
+
//________________________________________________________________________
Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
{
-
- if (centbin==0){
+ // Get phi mean.
+
+ if (centbin==0) {
Double_t PhiMean[9]={0.036,
0.021,
0.0121,
0.0022,
0.001};
return PhiMean[pbin];
- }else if(centbin==1){
+ } else if (centbin==1) {
Double_t PhiMean[9]={0.036,
0.021,
0.0121,
0.0022,
0.001};
return PhiMean[pbin];
- }else if(centbin==2){
+ } else if (centbin==2) {
Double_t PhiMean[9]={0.036,
0.021,
0.0121,
0.0022,
0.001};
return PhiMean[pbin];
- }else if(centbin==3){
+ } else if (centbin==3) {
Double_t PhiMean[9]={0.036,
0.021,
0.0121,
0.0022,
0.001};
return PhiMean[pbin];
- }else if(centbin==4){
+ } else if (centbin==4) {
Double_t PhiMean[9]={0.036,
0.021,
0.0127,
0.0028,
0.0018};
return PhiMean[pbin]*(-1.);
- }else if(centbin==5){
+ } else if (centbin==5) {
Double_t PhiMean[9]={0.036,
0.021,
0.0127,
0.0028,
0.0018};
return PhiMean[pbin]*(-1.);
- }else if(centbin==6){
+ } else if (centbin==6) {
Double_t PhiMean[9]={0.036,
0.021,
0.0127,
0.0028,
0.0018};
return PhiMean[pbin]*(-1.);
- }else if(centbin==7){
+ } else if (centbin==7) {
Double_t PhiMean[9]={0.036,
0.021,
0.0127,
}
return 0;
-
}
+
//________________________________________________________________________
Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
{
+ // Get phi sigma.
- if (centbin==0){
+ if (centbin==0) {
Double_t PhiSigma[9]={0.0221,
0.0128,
0.0074,
0.0049,
0.0045};
return 2.*PhiSigma[pbin];
- }else if(centbin==1){
+ } else if (centbin==1) {
Double_t PhiSigma[9]={0.0217,
0.0120,
0.0076,
0.0054,
0.0045};
return 2.*PhiSigma[pbin];
- }else if(centbin==2){
+ } else if (centbin==2) {
Double_t PhiSigma[9]={0.0211,
0.0124,
0.0080,
0.0054,
0.0047};
return 2.*PhiSigma[pbin];
- }else if(centbin==3){
+ } else if (centbin==3) {
Double_t PhiSigma[9]={0.0215,
0.0124,
0.0082,
0.0055,
0.0047};
return 2.*PhiSigma[pbin];
- }else if(centbin==4){
+ } else if (centbin==4) {
Double_t PhiSigma[9]={0.0199,
0.0108,
0.0072,
0.0049,
0.0045};
return 2.*PhiSigma[pbin];
- }else if(centbin==5){
+ } else if (centbin==5) {
Double_t PhiSigma[9]={0.0200,
0.0110,
0.0074,
0.0052,
0.0045};
return 2.*PhiSigma[pbin];
- }else if(centbin==6){
+ } else if (centbin==6) {
Double_t PhiSigma[9]={0.0202,
0.0113,
0.0077,
0.0055,
0.0050};
return 2.*PhiSigma[pbin];
- }else if(centbin==7){
+ } else if (centbin==7) {
Double_t PhiSigma[9]={0.0205,
0.0113,
0.0080,
}
return 0;
-
}
//________________________________________________________________________
AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
if (!handler) {
- AliError("Input handler not available!");
+ AliError(Form("%s: Input handler not available!", GetName()));
return;
}
} else if (handler->InheritsFrom("AliAODInputHandler")) {
fOutClusters = new TClonesArray("AliAODCaloCluster");
} else {
- AliError("Input handler not recognized!");
+ AliError(Form("%s: Input handler not recognized!", GetName()));
return;
}
fOutClusters->SetName(fOutCaloName);
+ if (!fCreateHisto)
+ return;
+
OpenFile(1);
- fOutputList = new TList();
- fOutputList->SetOwner();
+ fOutput = new TList();
+ fOutput->SetOwner();
- if (!fCreateHisto) {
- PostData(1, fOutputList);
- return;
- }
+ TString name;
for(Int_t icent=0; icent<8; ++icent) {
for(Int_t ipt=0; ipt<9; ++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]);
+ for(Int_t ieta=0; ieta<2; ++ieta) {
+ name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
+ fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, 200, -0.1, 0.1, 400, -0.2, 0.2);
+ fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
+ }
}
- if(icent<4){
- TString name(Form("fHistMatchEvsP_%i",icent));
+ for(Int_t itrk=0; itrk<4; ++itrk) {
+ name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
+ fHistNCellsEnergy[icent][itrk] = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
+ fOutput->Add(fHistNCellsEnergy[icent][itrk]);
+ }
+
+ name = Form("fHistEsubPch_%i",icent);
+ fHistEsubPch[icent]=new TH1F(name, name, 400, 0., 100.);
+ fOutput->Add(fHistEsubPch[icent]);
+
+ name = Form("fHistEsubPchRat_%i",icent);
+ fHistEsubPchRat[icent]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
+ fOutput->Add(fHistEsubPchRat[icent]);
+
+ if (icent<4) {
+ name = Form("fHistMatchEvsP_%i",icent);
fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
- fOutputList->Add(fHistMatchEvsP[icent]);
+ fOutput->Add(fHistMatchEvsP[icent]);
name = Form("fHistMatchdRvsEP_%i",icent);
fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
- fOutputList->Add(fHistMatchdRvsEP[icent]);
+ fOutput->Add(fHistMatchdRvsEP[icent]);
- for(Int_t itrk=0; itrk<3; ++itrk){
- name = Form("fHistEsubPch_%i_%i",icent,itrk);
- fHistEsubPch[icent][itrk]=new TH1F(name, name, 400, 0., 100.);
- fHistEsubPch[icent][itrk]->Sumw2();
- fOutputList->Add(fHistEsubPch[icent][itrk]);
-
- name = Form("fHistEsubPchRat_%i_%i",icent,itrk);
- fHistEsubPchRat[icent][itrk]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
- fOutputList->Add(fHistEsubPchRat[icent][itrk]);
- }
+ name = Form("fHistNMatchEnergy_%i",icent);
+ fHistNMatchEnergy[icent] = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
+ fOutput->Add(fHistNMatchEnergy[icent]);
}
}
fHistCentrality = new TH1F("fHistCentrality", "Centrality", 100, 0, 100);
-
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);
- fHistNMatchCent_trk = new TH2F("NMatchesCent_trk", "NMatchesCent_trk", 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);
- fOutputList->Add(fHistNMatchCent_trk);
- fOutputList->Add(fHistCentrality);
-
- PostData(1, fOutputList);
+ fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent", 100, 0, 100, 11, -0.5, 10.5);
+
+ fOutput->Add(fHistNclusMatchvsCent);
+ fOutput->Add(fHistNclusvsCent);
+ fOutput->Add(fHistEbefore);
+ fOutput->Add(fHistEafter);
+ fOutput->Add(fHistEoPCent);
+ fOutput->Add(fHistNMatchCent);
+ fOutput->Add(fHistCentrality);
+
+ PostData(1, fOutput);
}
//________________________________________________________________________
-void AliHadCorrTask::UserExec(Option_t *)
+void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
{
- // Execute per event.
+ // Do the loop over matched tracks for cluster emccluster.
+
+ AliVCluster *cluster = emccluster->GetCluster();
+ Int_t iClus = emccluster->IdInCollection();
+ Double_t energyclus = cluster->E();
+
+ // loop over matched tracks
+ const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
+ for (Int_t i = 0; i < Ntrks; ++i) {
+ Int_t iTrack = emccluster->GetMatchedObjId(i);
+ Double_t dR = emccluster->GetMatchedObjDistance(i);
+
+ AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
+ if (!emctrack)
+ continue;
+ AliVTrack *track = emctrack->GetTrack();
+ if (!track)
+ continue;
+ if (!AcceptTrack(track))
+ continue;
+
+ // check if track also points to cluster
+ if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
+ continue;
+
+ Double_t etadiff = 999;
+ Double_t phidiff = 999;
+ AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
+
+ Double_t mom = track->P();
+ Int_t mombin = GetMomBin(mom);
+ Int_t centbinch = fCentBin;
+ if (track->Charge()<0)
+ centbinch += 4;
+
+ if (fCreateHisto) {
+ Int_t etabin = 0;
+ if(track->Eta() > 0) etabin=1;
+ fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
+ }
+
+ Double_t etaCut = 0.0;
+ Double_t phiCutlo = 0.0;
+ Double_t phiCuthi = 0.0;
+
+ if (fPhiMatch > 0) {
+ phiCutlo = -fPhiMatch;
+ phiCuthi = +fPhiMatch;
+ } else {
+ phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
+ phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
+ }
+
+ if (fEtaMatch > 0) {
+ etaCut = fEtaMatch;
+ } else {
+ etaCut = GetEtaSigma(mombin);
+ }
+
+ if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
+ ++Nmatches;
+ totalTrkP += track->P();
+
+ if (fCreateHisto) {
+ if (fHadCorr > 1 && mombin > -1) {
+ fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
+ }
+ }
+ }
+ }
+}
+
+//________________________________________________________________________
+Bool_t AliHadCorrTask::Run()
+{
+ // Run the hadronic correction
// post output in event if not yet present
if (!(InputEvent()->FindListObject(fOutCaloName)))
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");
- am->LoadBranch("Centrality");
-
- TList *l = InputEvent()->GetList();
-
- // 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;
+ if (esdMode) { // optimization in case autobranch loading is off
+ AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
+ if (fCaloName == "CaloClusters")
+ am->LoadBranch("CaloClusters");
+ if (fTracksName == "Tracks")
+ am->LoadBranch("Tracks");
+ am->LoadBranch("Centrality.");
}
-
- Int_t centbin = GetCentBin(cent);
-
- if (fCreateHisto)
- fHistCentrality->Fill(cent);
- // get input collections
- TClonesArray *tracks = 0;
- TClonesArray *clus = 0;
-
- 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();
-
- if (fDoTrackClus) {
- for(Int_t t = 0; t<Ntrks; ++t) {
- AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(t));
- if (!track)
- continue;
- Int_t Nmatches = 0;
- Double_t dEtaMin = 1e9;
- Double_t dPhiMin = 1e9;
- Int_t imin = -1;
- for(Int_t i=0; i < Nclus; ++i) {
- AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(i));
- if (!c)
- continue;
- Double_t etadiff=999;
- Double_t phidiff=999;
- AliPicoTrack::GetEtaPhiDiff(track,c,phidiff,etadiff);
- Double_t dR = TMath::Sqrt(etadiff*etadiff+phidiff*phidiff);
- Double_t dRmin = TMath::Sqrt(dEtaMin*dEtaMin+dPhiMin*dPhiMin);
- if(dR > 25)
- continue;
- if(dR<dRmin){
- dEtaMin = etadiff;
- dPhiMin = phidiff;
- imin = i;
- }
- if (TMath::Abs(phidiff)<0.05 && TMath::Abs(etadiff)<0.025) { // pp cuts!!!
- Nmatches++;
- }
- }
-
- if (fCreateHisto)
- fHistNMatchCent_trk->Fill(cent,Nmatches);
-
- track->SetEMCALcluster(imin);
- }
-
+ if (fCreateHisto) {
+ fHistCentrality->Fill(fCent);
}
+ // loop over all clusters
+ const Int_t Nclus = fCaloClusters->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())
+ AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
+ if (!emccluster)
continue;
- // make primary particle out of cluster
- TLorentzVector nPart;
- c->GetMomentum(nPart, vertex);
-
- Double_t etclus = nPart.Pt();
- if (etclus<fMinPt)
+ AliVCluster *cluster = emccluster->GetCluster();
+ if (!cluster)
+ continue;
+ if (!AcceptCluster(cluster))
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;
- }
-
- Double_t mom = track->P();
- Int_t mombin = GetMomBin(mom);
- Int_t centbinch = centbin;
- if (track->Charge() == -1 || track->Charge() == 255)
- centbinch += 4;
-
- if (fCreateHisto) {
- if (fHadCorr > 1 && mombin > -1) {
- fHistMatchEtaPhi[centbinch][mombin]->Fill(etadiff, phidiff);
- fHistMatchdRvsEP[centbin]->Fill(dR, energyclus / mom);
- }
- }
-
- Double_t EtaCut = 0.0;
- Double_t PhiCutlo = 0.0;
- Double_t PhiCuthi = 0.0;
- if(fPhiMatch > 0){
- PhiCutlo = -fPhiMatch;
- PhiCuthi = fPhiMatch;
- }
- else {
- PhiCutlo = GetPhiMean(mombin,centbinch)-GetPhiSigma(mombin,centbin);
- PhiCuthi = GetPhiMean(mombin,centbinch)+GetPhiSigma(mombin,centbin);
- }
-
- if(fEtaMatch > 0){
- EtaCut = fEtaMatch;
- }
- else {
- EtaCut = GetEtaSigma(mombin);
- }
-
- if ((phidiff < PhiCuthi && phidiff > PhiCutlo) && TMath::Abs(etadiff) < EtaCut) {
- if((fDoTrackClus && (track->GetEMCALcluster()) == iClus) || !fDoTrackClus){
- ++Nmatches;
- totalTrkP += track->P();
- }
- }
- }
-
- // store closest track
- c->SetEmcCpvDistance(imin);
- c->SetTrackDistance(dPhiMin, dEtaMin);
-
- if(fCreateHisto) {
- fHistNclusvsCent->Fill(cent);
- fHistEbefore->Fill(cent, energyclus);
- fHistNMatchCent->Fill(cent, Nmatches);
-
- if(Nmatches > 0)
- fHistNclusMatchvsCent->Fill(cent);
- }
+ Double_t energyclus = 0;
+ if (fCreateHisto)
+ fHistEbefore->Fill(fCent, cluster->E());
// 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;
- Double_t Esub = (fHadCorr-1)*totalTrkP;
- if (Esub > energyclus)
- Esub = energyclus;
-
- if(fCreateHisto) {
- fHistEoPCent->Fill(cent, EoP);
- fHistMatchEvsP[centbin]->Fill(energyclus, EoP);
-
- if (Nmatches == 1) fHistEsubPchRat[centbin][0]->Fill(totalTrkP,Esub/totalTrkP);
- else if(Nmatches == 2) fHistEsubPchRat[centbin][1]->Fill(totalTrkP,Esub/totalTrkP);
- else fHistEsubPchRat[centbin][2]->Fill(totalTrkP,Esub/totalTrkP);
-
- if(Nmatches==1) fHistEsubPch[centbin][0]->Fill(totalTrkP,Esub);
- else if(Nmatches==2) fHistEsubPch[centbin][1]->Fill(totalTrkP,Esub);
- else fHistEsubPch[centbin][2]->Fill(totalTrkP,Esub);
- }
- }
- 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);
- Int_t centbinch = centbin;
- if (t->Charge() == -1 || t->Charge() == 255)
- centbinch += 4;
-
- if (fCreateHisto) {
- fHistMatchEtaPhi[centbinch][mombin]->Fill(dEtaMin, dPhiMin);
-
- if (mom > 0) {
- fHistMatchEvsP[centbin]->Fill(energyclus, energyclus / mom);
- fHistEoPCent->Fill(cent, energyclus / mom);
- fHistMatchdRvsEP[centbin]->Fill(dRmin, energyclus / mom);
- }
- }
-
- Double_t EtaCut = 0.0;
- Double_t PhiCutlo = 0.0;
- Double_t PhiCuthi = 0.0;
- if(fPhiMatch > 0) {
- PhiCutlo = -fPhiMatch;
- PhiCuthi = fPhiMatch;
- }
- else {
- PhiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, centbin);
- PhiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, centbin);
- }
-
- if(fEtaMatch > 0) {
- EtaCut = fEtaMatch;
- }
- else {
- EtaCut = GetEtaSigma(mombin);
- }
-
- if ((dPhiMin<PhiCuthi && dPhiMin>PhiCutlo) && TMath::Abs(dEtaMin) < EtaCut) {
-
- if((fDoTrackClus && (t->GetEMCALcluster()) == iClus) || !fDoTrackClus){
- energyclus -= fHadCorr * t->P();
- }
- }
- }
- }
+ if (fHadCorr > 1)
+ energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);
+ else
+ energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);
}
if (energyclus < 0)
energyclus = 0;
- if (fCreateHisto)
- 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)));
+ AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
+ if (!ec)
+ continue;
+ oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
+ } else {
+ AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
+ if (!ac)
+ continue;
+ oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
}
oc->SetE(energyclus);
+
++clusCount;
+
+ if (fCreateHisto)
+ fHistEafter->Fill(fCent, energyclus);
}
}
+
+ return kTRUE;
+}
- if (fCreateHisto)
- PostData(1, fOutputList);
+//________________________________________________________________________
+Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr)
+{
+ // Apply the hadronic correction with one track only.
+
+ AliVCluster *cluster = emccluster->GetCluster();
+ Double_t energyclus = cluster->E();
+ Int_t iMin = emccluster->GetMatchedObjId();
+ if (iMin < 0)
+ return energyclus;
+
+ AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
+ if (!emctrack)
+ return energyclus;
+
+ // check if track also points to cluster
+ Int_t cid = emctrack->GetMatchedObjId();
+ if (fDoTrackClus && (cid!=emccluster->IdInCollection()))
+ return energyclus;
+
+ AliVTrack *track = emctrack->GetTrack();
+ if (!track)
+ return energyclus;
+
+ Double_t mom = track->P();
+ if (mom == 0)
+ return energyclus;
+
+ Double_t dRmin = emccluster->GetMatchedObjDistance();
+ Double_t dEtaMin = 1e9;
+ Double_t dPhiMin = 1e9;
+ AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
+
+ Int_t mombin = GetMomBin(mom);
+ Int_t centbinch = fCentBin;
+ if (track->Charge()<0)
+ centbinch += 4;
+
+ // plot some histograms if switched on
+ if (fCreateHisto) {
+ Int_t etabin = 0;
+ if(track->Eta() > 0)
+ etabin = 1;
+
+ fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
+
+ if (mom > 0) {
+ fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
+ fHistEoPCent->Fill(fCent, energyclus / mom);
+ fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
+ }
+ }
+
+ // define eta/phi cuts
+ Double_t etaCut = 0.0;
+ Double_t phiCutlo = 0.0;
+ Double_t phiCuthi = 0.0;
+ if (fPhiMatch > 0) {
+ phiCutlo = -fPhiMatch;
+ phiCuthi = +fPhiMatch;
+ }
+ else {
+ phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
+ phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
+ }
+ if(fEtaMatch > 0) {
+ etaCut = fEtaMatch;
+ }
+ else {
+ etaCut = GetEtaSigma(mombin);
+ }
+
+ // apply the correction if the track is in the eta/phi window
+ if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
+ energyclus -= hadCorr * mom;
+ }
+
+ return energyclus;
}
//________________________________________________________________________
-void AliHadCorrTask::Terminate(Option_t *)
+Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
{
- // Nothing to be done.
+ // Apply the hadronic correction with all tracks.
+
+ AliVCluster *cluster = emccluster->GetCluster();
+
+ Double_t energyclus = cluster->E();
+ Double_t cNcells = cluster->GetNCells();
+
+ Double_t totalTrkP = 0.0; // count total track momentum
+ Int_t Nmatches = 0; // count total number of matches
+
+ // do the loop over the matched tracks and get the number of matches and the total momentum
+ DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
+
+ if(fCreateHisto && Nmatches == 0) {
+ fHistNclusvsCent->Fill(fCent);
+ fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
+ }
+
+ if (totalTrkP <= 0)
+ return energyclus;
+
+ Double_t Esub = hadCorr * totalTrkP;
+
+ Double_t clusEexcl = fEexclCell * cNcells;
+
+ if (energyclus < clusEexcl)
+ clusEexcl = energyclus;
+
+ if (Esub > energyclus)
+ Esub = energyclus;
+
+ // applying Peter's proposed algorithm
+ // never subtract the full energy of the cluster
+ if ((energyclus - Esub) < clusEexcl)
+ Esub = (energyclus - clusEexcl);
+
+ Double_t EoP = energyclus / totalTrkP;
+
+ // plot some histograms if switched on
+ if (fCreateHisto) {
+ fHistNMatchCent->Fill(fCent, Nmatches);
+ fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
+
+ if(Nmatches > 0)
+ fHistNclusMatchvsCent->Fill(fCent);
+
+ if(Nmatches == 1)
+ fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);
+ else if(Nmatches == 2)
+ fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
+ else if(Nmatches > 2)
+ fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
+
+ if (EoP > 0) {
+ fHistEoPCent->Fill(fCent, EoP);
+ fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
+ }
+
+ if (Nmatches == 1) {
+ Int_t iMin = emccluster->GetMatchedObjId();
+ AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
+ AliVTrack *track = emctrack->GetTrack();
+ Int_t centbinchm = fCentBin;
+ if (track->Charge()<0)
+ centbinchm += 4;
+
+ fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
+ fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
+ }
+ }
+
+ // apply the correction
+ energyclus -= Esub;
+
+ return energyclus;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff