22 #include <TClonesArray.h>
26 #include <THnSparse.h>
28 #include <TLorentzVector.h>
29 #include <TParameter.h>
30 #include <TParticle.h>
33 #include <TObjArray.h>
35 #include <AliVCluster.h>
36 #include <AliVParticle.h>
39 #include "AliAnalysisTask.h"
40 #include "AliAnalysisManager.h"
42 #include "AliESDEvent.h"
43 #include "AliESDInputHandler.h"
44 #include "AliESDtrackCuts.h"
45 #include "AliAODEvent.h"
46 #include "AliAODHandler.h"
56 #include "AliESDpid.h"
57 #include "AliAODPid.h"
58 #include "AliPIDResponse.h"
59 #include "AliKFParticle.h"
60 #include "AliKFVertex.h"
96 fHistClusEovPnonlin(),
97 fHistClusEovPHadCorr()
127 fHistClusEovPnonlin(0),
128 fHistClusEovPHadCorr(0)
151 AliVEventHandler *inputHandler =
dynamic_cast<AliVEventHandler *
>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
152 if(!TString(inputHandler->IsA()->GetName()).CompareTo(
"AliAODInputHandler")){
157 printf(
"Analysis Mode: %s Analysis\n",
IsAODanalysis() ?
"AOD" :
"ESD");
167 while ((obj = next())) {
201 while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
202 groupname = clusCont->GetName();
205 fHistClusEovP =
new TH1F(
"ClusEovP",
"Cluster EovP",130,0,1.3);
208 fInvmassULS =
new TH1F(
"fInvmassULS",
"Inv mass of ULS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
209 fInvmassLS =
new TH1F(
"fInvmassLS",
"Inv mass of LS (e,e) for pt^{e}>1; mass(GeV/c^2); counts;", 1000,0,1.0);
211 for (Int_t cent = 0; cent <
fNcentBins; cent++) {
212 histname = TString::Format(
"%s/histClusterEnergy_%d", groupname.Data(), cent);
213 histtitle = TString::Format(
"%s;#it{E}_{cluster} (GeV);counts", histname.Data());
216 histname = TString::Format(
"%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
217 histtitle = TString::Format(
"%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
220 histname = TString::Format(
"%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
221 histtitle = TString::Format(
"%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
224 histname = TString::Format(
"%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
225 histtitle = TString::Format(
"%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
228 histname = TString::Format(
"%s/histClusterPhi_%d", groupname.Data(), cent);
229 histtitle = TString::Format(
"%s;#it{#phi}_{custer};counts", histname.Data());
232 histname = TString::Format(
"%s/histClusterEta_%d", groupname.Data(), cent);
233 histtitle = TString::Format(
"%s;#it{#eta}_{custer};counts", histname.Data());
236 histname = TString::Format(
"%s/histNClusters_%d", groupname.Data(), cent);
237 histtitle = TString::Format(
"%s;number of clusters;events", histname.Data());
260 for (Int_t cent = 0; cent <
fNcentBins; cent++) {
261 histname = TString::Format(
"%s/histCellEnergyvsAbsId_%d", groupname.Data(), cent);
262 histtitle = TString::Format(
"%s;cell abs. ID;#it{E}_{cell} (GeV);counts", histname.Data());
265 histname = TString::Format(
"%s/histNCells_%d", groupname.Data(), cent);
266 histtitle = TString::Format(
"%s;number of cells;events", histname.Data());
289 while ((partCont = static_cast<AliParticleContainer*>(next()))) {
290 groupname = partCont->GetName();
292 for (Int_t cent = 0; cent <
fNcentBins; cent++) {
293 histname = TString::Format(
"%s/histTrackPt_%d", groupname.Data(), cent);
294 histtitle = TString::Format(
"%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
297 histname = TString::Format(
"%s/histTrackPhi_%d", groupname.Data(), cent);
298 histtitle = TString::Format(
"%s;#it{#phi}_{track};counts", histname.Data());
301 histname = TString::Format(
"%s/histTrackEta_%d", groupname.Data(), cent);
302 histtitle = TString::Format(
"%s;#it{#eta}_{track};counts", histname.Data());
305 if (TClass(partCont->GetClassName()).InheritsFrom(
"AliVTrack")) {
306 histname = TString::Format(
"%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
307 histtitle = TString::Format(
"%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
310 histname = TString::Format(
"%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
311 histtitle = TString::Format(
"%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
314 histname = TString::Format(
"%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
315 histtitle = TString::Format(
"%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{p}_{T,track}^{vertex} - #it{p}_{T,track}^{EMCal} (GeV/#it{c});counts", histname.Data());
318 histname = TString::Format(
"%s/fHistEoverPvsP_%d", groupname.Data(), cent);
319 histtitle = TString::Format(
"%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
323 histname = TString::Format(
"%s/histNTracks_%d", groupname.Data(), cent);
324 histtitle = TString::Format(
"%s;number of tracks;events", histname.Data());
347 while ((jetCont = static_cast<AliJetContainer*>(next()))) {
348 groupname = jetCont->GetName();
351 fHistJetEovP =
new TH1F(
"JetEovP",
"HFE Jet EovP",200,0,5);
352 fHistJetEovPvPt =
new TH2F(
"JetEovPvPt",
"HFE Jet EovP vs Pt",100,0,1,100,0,100);
354 for (Int_t cent = 0; cent <
fNcentBins; cent++) {
355 histname = TString::Format(
"%s/histJetPt_%d", groupname.Data(), cent);
356 histtitle = TString::Format(
"%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
359 histname = TString::Format(
"%s/histJetArea_%d", groupname.Data(), cent);
360 histtitle = TString::Format(
"%s;#it{A}_{jet};counts", histname.Data());
363 histname = TString::Format(
"%s/histJetPhi_%d", groupname.Data(), cent);
364 histtitle = TString::Format(
"%s;#it{#phi}_{jet};counts", histname.Data());
367 histname = TString::Format(
"%s/histJetEta_%d", groupname.Data(), cent);
368 histtitle = TString::Format(
"%s;#it{#eta}_{jet};counts", histname.Data());
371 histname = TString::Format(
"%s/histNJets_%d", groupname.Data(), cent);
372 histtitle = TString::Format(
"%s;number of jets;events", histname.Data());
381 histname = TString::Format(
"%s/histJetCorrPt_%d", groupname.Data(), cent);
382 histtitle = TString::Format(
"%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
416 while ((jetCont = static_cast<AliJetContainer*>(next()))) {
417 groupname = jetCont->GetName();
422 for(
auto jet : jetCont->
accepted()) {
426 histname = TString::Format(
"%s/histJetPt_%d", groupname.Data(),
fCentBin);
429 histname = TString::Format(
"%s/histJetArea_%d", groupname.Data(),
fCentBin);
432 histname = TString::Format(
"%s/histJetPhi_%d", groupname.Data(),
fCentBin);
435 histname = TString::Format(
"%s/histJetEta_%d", groupname.Data(),
fCentBin);
439 histname = TString::Format(
"%s/histJetCorrPt_%d", groupname.Data(),
fCentBin);
442 histname = TString::Format(
"%s/histNJets_%d", groupname.Data(),
fCentBin);
448 Float_t corrPt = jet->Pt() - rhoVal * jet->Area();
449 TLorentzVector leadPart;
456 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
457 AliVParticle *track = jet->TrackAt(it, tracks->GetArray());
460 Double_t JetTrackP = track->P();
461 Double_t JetTrackPt = track->Pt();
464 Double_t dphi = TVector2::Phi_0_2pi(track->Phi() - jet->Phi());
465 Double_t deta = track->Eta() - jet->Eta();
466 Double_t dist = TMath::Sqrt(deta * deta + dphi * dphi);
469 AliVTrack *Vtrack =
dynamic_cast<AliVTrack*
>(track);
470 Int_t EMCalIndex = -1;
471 EMCalIndex = Vtrack->GetEMCALcluster();
472 if(EMCalIndex < 0)
continue;
474 AliVCluster *clustMatch = (AliVCluster*)
fVevent->GetCaloCluster(EMCalIndex);
475 if(!(clustMatch && clustMatch->IsEMCAL()))
continue;
477 Double_t clustMatchE = clustMatch->E();
478 Double_t clustMatchNonLinE = clustMatch->GetNonLinCorrEnergy();
479 Double_t clustMatchHadCorr = clustMatch->GetHadCorrEnergy();
504 while ((partCont = static_cast<AliParticleContainer*>(next()))) {
505 groupname = partCont->GetName();
507 for(
auto part : partCont->
accepted()) {
511 histname = TString::Format(
"%s/histTrackPt_%d", groupname.Data(),
fCentBin);
514 histname = TString::Format(
"%s/histTrackPhi_%d", groupname.Data(),
fCentBin);
517 histname = TString::Format(
"%s/histTrackEta_%d", groupname.Data(),
fCentBin);
520 if (partCont->GetLoadedClass()->InheritsFrom(
"AliVTrack")) {
521 const AliVTrack* track =
static_cast<const AliVTrack*
>(part);
523 histname = TString::Format(
"%s/fHistDeltaEtaPt_%d", groupname.Data(),
fCentBin);
526 histname = TString::Format(
"%s/fHistDeltaPhiPt_%d", groupname.Data(),
fCentBin);
529 histname = TString::Format(
"%s/fHistDeltaPtvsPt_%d", groupname.Data(),
fCentBin);
532 Double_t TrackP = track->P();
538 Bool_t fFlagPhotonicElec = kFALSE;
544 Int_t iCluster = track->GetEMCALcluster();
549 Double_t ClusterE = cluster->E();
550 Double_t ClusterNonLinE = cluster->GetNonLinCorrEnergy();
551 Double_t ClusterHadCorr = cluster->GetHadCorrEnergy();
553 Double_t
sigma = 0.0;
556 sigma =
fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
562 histname = TString::Format(
"%s/fHistEoverPvsP_%d", groupname.Data(),
fCentBin);
570 histname = TString::Format(
"%s/histNTracks_%d", groupname.Data(),
fCentBin);
585 while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
586 groupname = clusCont->GetName();
588 for(
auto cluster : clusCont->
all()) {
589 if (!cluster)
continue;
591 if (cluster->GetIsExotic()) {
592 histname = TString::Format(
"%s/histClusterEnergyExotic_%d", groupname.Data(),
fCentBin);
598 for(
auto cluster : clusCont->
accepted()) {
599 if (!cluster)
continue;
603 cluster->GetMomentum(nPart,
fVertex);
605 histname = TString::Format(
"%s/histClusterEnergy_%d", groupname.Data(),
fCentBin);
608 histname = TString::Format(
"%s/histClusterNonLinCorrEnergy_%d", groupname.Data(),
fCentBin);
611 histname = TString::Format(
"%s/histClusterHadCorrEnergy_%d", groupname.Data(),
fCentBin);
614 histname = TString::Format(
"%s/histClusterPhi_%d", groupname.Data(),
fCentBin);
617 histname = TString::Format(
"%s/histClusterEta_%d", groupname.Data(),
fCentBin);
621 histname = TString::Format(
"%s/histNClusters_%d", groupname.Data(),
fCentBin);
636 const Short_t ncells =
fCaloCells->GetNumberOfCells();
642 for (Short_t pos = 0; pos < ncells; pos++) {
644 Short_t absId =
fCaloCells->GetCellNumber(pos);
668 UInt_t evSelMask=((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
670 fVevent =
dynamic_cast<AliVEvent*
>(InputEvent());
672 printf(
"ERROR: fVEvent not available\n");
676 fESD =
dynamic_cast<AliESDEvent*
>(InputEvent());
682 fAOD =
dynamic_cast<AliAODEvent*
>(InputEvent());
688 const AliVVertex *pVtx =
fVevent->GetPrimaryVertex();
692 printf(
"ERROR: fpidResponse not available\n");
697 ntracks =
fVevent->GetNumberOfTracks();
700 Double_t Zvertex = -100, Xvertex = -100, Yvertex = -100;
701 Double_t NcontV = pVtx->GetNContributors();
711 Int_t numberofvertices = 100;
712 if(
fAOD) numberofvertices =
fAOD->GetNumberOfVertices();
713 Double_t listofmotherkink[numberofvertices];
714 Int_t numberofmotherkink = 0;
717 for(Int_t ivertex=0; ivertex < numberofvertices; ivertex++) {
718 AliAODVertex *aodvertex =
fAOD->GetVertex(ivertex);
719 if(!aodvertex)
continue;
720 if(aodvertex->GetType()==AliAODVertex::kKink) {
721 AliAODTrack *mother = (AliAODTrack *) aodvertex->GetParent();
722 if(!mother)
continue;
723 Int_t idmother = mother->GetID();
724 listofmotherkink[numberofmotherkink] = idmother;
725 numberofmotherkink++;
747 Double_t clustE=-999, emcphi = -999, emceta=-999;
748 Nclust =
fVevent->GetNumberOfCaloClusters();
749 for(Int_t icl=0; icl<Nclust; icl++)
751 AliVCluster *clust = 0x0;
752 clust =
fVevent->GetCaloCluster(icl);
753 if(!clust) printf(
"ERROR: Could not receive cluster matched calibrated from track %d\n", icl);
755 if(clust && clust->IsEMCAL())
758 clust->GetPosition(emcx);
759 clustpos.SetXYZ(emcx[0],emcx[1],emcx[2]);
760 emcphi = clustpos.Phi();
761 emceta = clustpos.Eta();
777 Bool_t flagPhotonicElec = kFALSE;
778 Double_t ptAsso=-999.,
nsigma=-999.0;
779 Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
781 for(Int_t jTracks = 0; jTracks<
fVevent->GetNumberOfTracks(); jTracks++){
782 if(jTracks==itrack)
continue;
784 AliVParticle* VtrackAsso =
fVevent->GetTrack(jTracks);
786 printf(
"ERROR: Could not receive track %d\n", jTracks);
790 AliVTrack *trackAsso =
dynamic_cast<AliVTrack*
>(VtrackAsso);
791 if(!trackAsso)
continue;
794 AliAODTrack *atrackAsso =
dynamic_cast<AliAODTrack*
>(VtrackAsso);
795 if(!atrackAsso)
continue;
796 if(!atrackAsso->TestFilterMask(AliAODTrack::kTrkTPCOnly))
continue;
797 if(atrackAsso->GetTPCNcls() < 70)
continue;
798 if(!(atrackAsso->GetStatus()&AliESDtrack::kTPCrefit))
continue;
799 if(!(atrackAsso->GetStatus()&AliESDtrack::kITSrefit))
continue;
803 ptAsso = trackAsso->Pt();
804 Int_t chargeAsso = trackAsso->Charge();
805 Int_t
charge = track->Charge();
807 if(ptAsso <0.2)
continue;
808 if(trackAsso->Eta()<-0.9 || trackAsso->Eta()>0.9)
continue;
809 if(nsigma < -3 || nsigma > 3)
continue;
811 Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
812 if(charge>0) fPDGe1 = -11;
813 if(chargeAsso>0) fPDGe2 = -11;
815 if(charge == chargeAsso) fFlagLS = kTRUE;
816 if(charge != chargeAsso) fFlagULS = kTRUE;
818 AliKFParticle::SetField(
fVevent->GetMagneticField());
820 AliKFParticle ge1 = AliKFParticle(*track, fPDGe1);
821 AliKFParticle ge2 = AliKFParticle(*trackAsso, fPDGe2);
822 AliKFParticle recg(ge1, ge2);
824 if(recg.GetNDF()<1)
continue;
825 Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
826 if(TMath::Sqrt(TMath::Abs(chi2recg))>3.)
continue;
828 Double_t
mass=-999., width = -999;
830 MassCorrect = recg.GetMass(mass,width);
832 if(fFlagLS && track->Pt()>1)
fInvmassLS->Fill(mass);
833 if(fFlagULS && track->Pt()>1)
fInvmassULS->Fill(mass);
835 if(mass<
fInvmassCut && fFlagULS && !flagPhotonicElec){
836 flagPhotonicElec = kTRUE;
839 fFlagPhotonicElec = flagPhotonicElec;
848 cout<<
"#######################"<<endl;
849 cout<<
"#### Task Finished ####"<<endl;
850 cout<<
"#######################"<<endl;
virtual ~AliAnalysisTaskEmcalJetHF()
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
const TString & GetRhoName() const
void AllocateClusterHistograms()
AliAODEvent * fAOD
ESD object.
void Terminate(Option_t *option)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Declaration of class AliTLorentzVector.
Double_t fMinBinPt
min pt in histograms
void AllocateTrackHistograms()
void AllocateCellHistograms()
Int_t fCentBin
!event centrality bin
TH1F * fHistClusEovPHadCorr
void SelectPhotonicElectron(Int_t itrack, AliVTrack *track, Bool_t &fFlagPhotonicElec)
TH2F * fM02EovP
M20 vs E/p.
THistManager fHistManager
Histogram manager.
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
particle/track collection array
const AliClusterIterableContainer all() const
THashList * CreateHistoGroup(const char *groupname, const char *parent="/")
AliAnalysisTaskEmcalJetHF()
AliParticleContainer * GetParticleContainer() const
void GetLeadingHadronMomentum(TLorentzVector &mom, const AliEmcalJet *jet) const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
THashList * GetListOfHistograms() const
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
Double_t Phi_0_2pi() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliVCluster * GetAcceptCluster(Int_t i) const
const AliClusterIterableContainer accepted() const
TString fCaloCellsName
name of calo cell collection
Bool_t IsAODanalysis() const
void AllocateJetHistograms()
Implementation of a HFE jet analysis task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
TObjArray fJetCollArray
jet collection array
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1F * fHistClusEovPnonlin
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t fVertex[3]
!event vertex
void UserCreateOutputObjects()
TH1F * fInvmassLS
M20 vs E/p.
TH1F * fEMCTrketa
Invmass of ULS pairs.
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
AliESDEvent * fESD
event object
const AliParticleIterableContainer accepted() const
TH1F * fHistJetEovP
EMC trk phi.
void UserCreateOutputObjects()
Double_t fEventCounter
pid response
const AliJetIterableContainer accepted() const
TH1F * fEMCTrkphi
EMC trk eta.
Container structure for EMCAL clusters.
TH1F * fInvmassULS
Invmass of LS pairs.
AliPIDResponse * fpidResponse
AOD object.
Container for jet within the EMCAL jet framework.
Int_t fNbins
no. of pt bins
Declaration of class AliAnalysisTaskEmcalJetHF.