//
// 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 "AliAnalysisManager.h"
-#include "AliAODCaloCluster.h"
#include "AliAODEvent.h"
-#include "AliESDEvent.h"
+#include "AliAODCaloCluster.h"
#include "AliESDCaloCluster.h"
-#include "AliCentrality.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),
fHistEoPCent(0),
fHistNMatchCent(0),
fHistNMatchCent_trk(0),
- fHistCentrality(0)
+ fHistCentrality(0),
+ fHistNoMatchEtaPhi(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"),
+ AliAnalysisTaskEmcal(name, kFALSE),
fOutCaloName("CaloClustersCorr"),
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),
fHistEoPCent(0),
fHistNMatchCent(0),
fHistNMatchCent_trk(0),
- fHistCentrality(0)
+ fHistCentrality(0),
+ fHistNoMatchEtaPhi(0)
{
// Standard constructor.
- for(Int_t i=0; i<8; i++) {
+ for(Int_t i=0; i<8; i++) {
+ fHistEsubPch[i] = 0;
+ fHistEsubPchRat[i] = 0;
+ for(Int_t j=0; j<3; j++) {
+ }
+
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;
- }
}
+
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.";
}
//________________________________________________________________________
AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
- AliAnalysisTaskSE(name),
- fTracksName("Tracks"),
- fCaloName("CaloClusters"),
+ AliAnalysisTaskEmcal(name, histo),
fOutCaloName("CaloClustersCorr"),
fPhiMatch(0.05),
fEtaMatch(0.025),
fDoTrackClus(0),
fHadCorr(0),
- fMinPt(0.15),
- fCreateHisto(histo),
+ fEexclCell(0),
fOutClusters(0),
- fOutputList(0),
fHistNclusvsCent(0),
fHistNclusMatchvsCent(0),
fHistEbefore(0),
fHistEoPCent(0),
fHistNMatchCent(0),
fHistNMatchCent_trk(0),
- fHistCentrality(0)
+ fHistCentrality(0),
+ fHistNoMatchEtaPhi(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.";
-
- if (fCreateHisto)
- 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
{
//________________________________________________________________________
Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
{
-
- 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){
Double_t PhiMean[9]={0.036,
0.021,
return 0;
}
+
//________________________________________________________________________
Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
{
-
if (centbin==0){
Double_t PhiSigma[9]={0.0221,
0.0128,
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]);
+ }
}
+
+ 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){
- TString name(Form("fHistMatchEvsP_%i",icent));
+ 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]);
+
+ name = Form("fHistNMatchEnergy_%i",icent);
+ fHistNMatchEnergy[icent] = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
+ fOutput->Add(fHistNMatchEnergy[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]);
- }
}
}
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);
+ fHistNMatchCent_trk = new TH2F("NMatchesCent_trk", "NMatchesCent_trk", 100, 0, 100, 11, -0.5, 10.5);
+
+ fHistNoMatchEtaPhi = new TH2F("NoMatchEtaPhi","NoMatchEtaPhi",200,-1.0,1.0,90,1.0,4.0);
+
+ fOutput->Add(fHistNclusMatchvsCent);
+ fOutput->Add(fHistNclusvsCent);
+ fOutput->Add(fHistEbefore);
+ fOutput->Add(fHistEafter);
+ fOutput->Add(fHistEoPCent);
+ fOutput->Add(fHistNMatchCent);
+ fOutput->Add(fHistNMatchCent_trk);
+ fOutput->Add(fHistCentrality);
+ fOutput->Add(fHistNoMatchEtaPhi);
+
+ PostData(1, fOutput);
}
-//_____________________________________________________
-TString AliHadCorrTask::GetBeamType()
+//________________________________________________________________________
+void AliHadCorrTask::DoTrackClusLoop()
{
- // Get beam type : pp-AA-pA
- // ESDs have it directly, AODs get it from hardcoded run number ranges
-
- AliVEvent *event = InputEvent();
+ const Int_t Ntrks = fTracks->GetEntries();
- if (!event) {
- AliError("Couldn't retrieve event!");
- return "";
- }
+ // loop over all tracks
+ for (Int_t t = 0; t < Ntrks; ++t) {
+ AliEmcalParticle *emctrack = dynamic_cast<AliEmcalParticle*>(fTracks->At(t));
+ if (!emctrack)
+ continue;
- TString beamType;
+ AliVTrack *track = emctrack->GetTrack();
+ if (!track)
+ continue;
+ if (!AcceptTrack(track))
+ continue;
+
+ Int_t Nclus = emctrack->GetNumberOfMatchedObj();
+ Int_t Nmatches = 0;
+
+ // loop over matched clusters
+ for (Int_t i = 0; i < Nclus; ++i) {
+ Int_t c = emctrack->GetMatchedObjId(i);
+ AliEmcalParticle *emccluster = dynamic_cast<AliEmcalParticle*>(fCaloClusters->At(c));
+ if (!emccluster)
+ continue;
+
+ AliVCluster *cluster = emccluster->GetCluster();
+ if (!cluster)
+ continue;
+
+ Double_t etadiff = 999;
+ Double_t phidiff = 999;
+ AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
+
+ if (TMath::Abs(phidiff) < 0.050 && TMath::Abs(etadiff) < 0.025) // pp cuts!!!
+ Nmatches++;
+ }
- AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
- if (esd) {
- const AliESDRun *run = esd->GetESDRun();
- beamType = run->GetBeamType();
+ fHistNMatchCent_trk->Fill(fCent, Nmatches);
+
+ if(Nmatches == 0 && track->Pt() > 2.0)
+ fHistNoMatchEtaPhi->Fill(track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal());
}
- else
- {
- Int_t runNumber = event->GetRunNumber();
- if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
- (runNumber >= 166529 && runNumber <= 170593)) // LHC11h
- {
- beamType = "A-A";
+}
+
+//________________________________________________________________________
+void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
+{
+ // 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
+ 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 = dynamic_cast<AliEmcalParticle*>(fTracks->At(iTrack));
+ if (!emctrack)
+ continue;
+
+ AliVTrack *track = emctrack->GetTrack();
+ if (!track)
+ continue;
+ if (!AcceptTrack(track))
+ 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() == -1 || track->Charge() == 255)
+ 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
- {
- beamType = "p-p";
+ else {
+ PhiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
+ PhiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
+ }
+
+ if (fEtaMatch > 0) {
+ EtaCut = fEtaMatch;
+ }
+ else {
+ EtaCut = GetEtaSigma(mombin);
}
- }
- return beamType;
+ if ((phidiff < PhiCuthi && phidiff > PhiCutlo) && TMath::Abs(etadiff) < EtaCut) {
+ if ((fDoTrackClus && (track->GetEMCALcluster()) == iClus) || !fDoTrackClus) {
+ ++Nmatches;
+ totalTrkP += track->P();
+
+ if (fCreateHisto) {
+ if (fHadCorr > 1 && mombin > -1) {
+ fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
+ }
+ }
+ }
+ }
+ }
}
//________________________________________________________________________
-void AliHadCorrTask::UserExec(Option_t *)
+Bool_t AliHadCorrTask::Run()
{
- // Execute per event.
+ // Run the hadronic correction
// post output in event if not yet present
if (!(InputEvent()->FindListObject(fOutCaloName)))
if (fTracksName == "Tracks")
am->LoadBranch("Tracks");
am->LoadBranch("Centrality");
-
- TList *l = InputEvent()->GetList();
-
- // get centrality
- Double_t cent = 99;
-
- if (GetBeamType() == "A-A") {
- AliCentrality *centrality = InputEvent()->GetCentrality();
-
- if (centrality)
- cent = centrality->GetCentralityPercentile("V0M");
- else
- cent = 99; // probably pp data
-
- if (cent < 0) {
- AliWarning(Form("Centrality negative: %f, assuming 99", cent));
- cent = 99;
- }
- }
-
- 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);
+ fHistCentrality->Fill(fCent);
- // loop over clusters and tracks
- const Int_t Nclus = clus->GetEntries();
- const Int_t Ntrks = tracks->GetEntries();
+ const Int_t Nclus = fCaloClusters->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 (fDoTrackClus && fCreateHisto)
+ DoTrackClusLoop();
+ // loop over all clusters
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 = dynamic_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;
+ 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)));
- }
+ if (esdMode)
+ oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*(dynamic_cast<AliESDCaloCluster*>(cluster)));
+ else
+ oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*(dynamic_cast<AliAODCaloCluster*>(cluster)));
+
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)
+{
+ AliVCluster *cluster = emccluster->GetCluster();
+
+ Double_t energyclus = cluster->E();
+ Int_t iMin = emccluster->GetMatchedObjId();
+
+ if (iMin < 0)
+ return energyclus;
+
+ AliEmcalParticle *emctrack = dynamic_cast<AliEmcalParticle*>(fTracks->At(iMin));
+ if (!emctrack)
+ 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() == -1 || track->Charge() == 255)
+ 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) {
+
+ if ((fDoTrackClus && (track->GetEMCALcluster()) == emccluster->IdInCollection()) || !fDoTrackClus){
+ energyclus -= hadCorr * mom;
+ }
+ }
+
+ return energyclus;
}
//________________________________________________________________________
-void AliHadCorrTask::Terminate(Option_t *)
+Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
{
- // Nothing to be done.
+ 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);
+
+ Double_t Esub = hadCorr * totalTrkP;
+
+ Double_t EoP = -1;
+ if (totalTrkP > 0)
+ EoP = energyclus / totalTrkP;
+
+ // Plot some histograms if switched on
+ if(fCreateHisto) {
+ fHistNclusvsCent->Fill(fCent);
+ fHistNMatchCent->Fill(fCent, Nmatches);
+ fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
+
+ if(Nmatches > 0)
+ fHistNclusMatchvsCent->Fill(fCent);
+
+ if(Nmatches == 0)
+ fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
+ else if(Nmatches == 1)
+ fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);
+ else if(Nmatches == 2)
+ fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
+ else
+ 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 = dynamic_cast<AliEmcalParticle*>(fTracks->At(iMin));
+ AliVTrack *track = emctrack->GetTrack();
+ Int_t centbinchm = fCentBin;
+ if (track->Charge() == -1 || track->Charge() == 255)
+ centbinchm += 4;
+
+ if (totalTrkP > 0)
+ fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
+
+ fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
+ }
+ }
+
+ if (totalTrkP <= 0)
+ return energyclus;
+
+ 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);
+
+ //apply the correction
+ energyclus -= Esub;
+
+ return energyclus;
}
+void AliHadCorrTask::Terminate(Option_t *)
+{
+ // Nothing to be done.
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff