From fae5f6bbc89f8301bedd98468909336967b211bc Mon Sep 17 00:00:00 2001 From: slindal Date: Fri, 8 Oct 2010 09:44:51 +0000 Subject: [PATCH] Updated gamma jet task. --- PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx | 314 ++++++++++++++------- PWG4/GammaConv/AliAnalysisTaskGammaJet.h | 54 +++- 2 files changed, 269 insertions(+), 99 deletions(-) diff --git a/PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx b/PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx index fa61bae7112..7a72a38f5e2 100644 --- a/PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx +++ b/PWG4/GammaConv/AliAnalysisTaskGammaJet.cxx @@ -1,7 +1,9 @@ +#include #include "TChain.h" #include "TTree.h" #include "TH1F.h" #include "TCanvas.h" +#include "TString.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" @@ -11,14 +13,12 @@ #include "AliESDCaloCluster.h" #include "AliESDInputHandler.h" - - +#include "AliAODPWG4ParticleCorrelation.h" #include "AliAODEvent.h" #include "AliAODHandler.h" #include "AliAODCaloCluster.h" #include "AliGammaConversionAODObject.h" - -#include +#include "AliAODJet.h" // Gamma - jet correlation analysis task // Authors: Svein Lindal @@ -36,35 +36,38 @@ AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet() : AliAnalysisTaskSE(), fHistPtEmcal(NULL), fHistPtJets(NULL), fHistGammaJets(NULL), - fDeltaAODFileName("") + fHistGammaJetsIso(NULL), + fMinPt(5.0), + fConeSize(0.3), + fPtThreshold(2.0), + fDeltaAODFileName(""), + fPhotons(NULL) { - // Constructor - - // Define input and output slots here - // Input slot #0 works with a TChain - //DefineInput(0, TChain::Class()); - // Output slot #0 id reserved by the base class for AOD - // Output slot #1 writes into a TH1 container - //DefineOutput(1, TList::Class()); + // Dummy Constructor } //________________________________________________________________________ -AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet(const char *name) : AliAnalysisTaskSE(name), +AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet(const char *name) : + AliAnalysisTaskSE(name), fOutputList(0), fHistPt(0), fHistPtPhos(0), fHistPtEmcal(0), fHistPtJets(0), fHistGammaJets(NULL), - fDeltaAODFileName("") + fHistGammaJetsIso(NULL), + fMinPt(0.0), + fConeSize(0.0), + fPtThreshold(0.0), + fDeltaAODFileName(""), + fPhotons(NULL) { // Constructor - // Define input and output slots here - // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); // Output slot #0 id reserved by the base class for AOD + // Output slot #1 writes into a TH1 container DefineOutput(1, TList::Class()); } @@ -72,39 +75,53 @@ AliAnalysisTaskGammaJet::AliAnalysisTaskGammaJet(const char *name) : AliAnalysis //________________________________________________________________________ void AliAnalysisTaskGammaJet::UserCreateOutputObjects() { - // Create histograms - // Called once - + //Create histograms add, to outputlist fOutputList = new TList(); fHistPt = new TH1F("fHistPt", "P_{T} distribution", 150, 0.1, 50); fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPt->SetMarkerStyle(kFullCircle); + fOutputList->Add(fHistPt); fHistPtPhos = new TH1F("fHistPtPhos", "P_{T} distribution", 150, 0.1, 50); fHistPtPhos->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPtPhos->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPtPhos->SetMarkerStyle(kFullCircle); + fOutputList->Add(fHistPtPhos); fHistPtEmcal = new TH1F("fHistPtEmcal", "P_{T} distribution", 150, 0.1, 50); fHistPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPtEmcal->SetMarkerStyle(kFullCircle); + fOutputList->Add(fHistPtEmcal); fHistPtJets = new TH1F("fHistPtJets", "P_{T} distribution", 150, 0.1, 50); fHistPtJets->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPtJets->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPtJets->SetMarkerStyle(kFullCircle); + fOutputList->Add(fHistPtJets); fHistGammaJets = new TH1F("fHistGammaJets", "fHistGammaJets", 200, -TMath::Pi(), 2*TMath::Pi()); - - fOutputList->Add(fHistPt); - fOutputList->Add(fHistPtPhos); - fOutputList->Add(fHistPtEmcal); - fOutputList->Add(fHistPtJets); fOutputList->Add(fHistGammaJets); + + fHistGammaJetsIso = new TH1F("fHistGammaJetsIso", "fHistGammaJetsIso", 200, -TMath::Pi(), 2*TMath::Pi()); + fOutputList->Add(fHistGammaJetsIso); + + //TNtuple * tuple = new TNtuple("fNtuple", "fNtuple", dPhi, + + + ///Create AOD branch + fPhotons = new TClonesArray("AliAODPWG4ParticleCorrelation", 0); + fPhotons->SetName("fPhotons"); + AddAODBranch("TClonesArray", &fPhotons); + + ///Isolation class + // fIsolation = new AliAnaParticleIsolation(); + // fIsolation->SetInputAODName("fPhotons"); + + } //________________________________________________________________________ @@ -112,114 +129,217 @@ void AliAnalysisTaskGammaJet::UserExec(Option_t *) { // Main loop // Called for each event + + //Clear stuff for new event + CleanUp(); - printf("in userexec \n"); - - // fESD = dynamic_cast(InputEvent()); - // if (!fESD) { - // printf("ERROR: fESD not available\n"); - // //return; - // AliAODEvent * fAOD = dynamic_cast(InputEvent()); - - // AliAODEvent * aodEvent = AODEvent(); - //AliESDEvent * esdEvent = dynamic_cast(InputEvent()); - - + ///Get AOD event AliAODEvent * aodEvent = GetAODEvent(); if(!aodEvent) { - cout << "No AOD event!!"<(aodEvent->FindListObject("GammaConv")); - if(!convGamma) - convGamma = GetConversionGammas(); - - if(!convGamma) { - printf("No convgamma"); - return; - } - - TClonesArray * jets = aodEvent->GetJets(); + //FillPWG4PartCorrBranch(convGamma, fPhotons, "ConvGamma"); + //fIsolation->MakeAnalysisFillAOD(); - for (Int_t iPhot = 0; iPhot < convGamma->GetEntriesFast(); iPhot++) { - AliGammaConversionAODObject * aodO = dynamic_cast(convGamma->At(iPhot)); - if (!aodO) { - printf("ERROR: Could not receive ga %d\n", iPhot); - continue; - } - - //if(aodO->E() < 0.2) continue; - - if (jets) { - for(int ij = 0; ij < jets->GetEntriesFast(); ij++) { - AliAODJet * jet = aodEvent->GetJet(ij); - if(jet) { - //cout << jet->E() << endl; - fHistPtJets->Fill(jet->E()); - cout << jet->Phi() << " " << aodO->Phi() << endl; - fHistGammaJets->Fill(jet->Phi() - aodO->Phi()); - - } - } - } - + ProcessConvGamma(aodEvent); + ProcessCalorimeters(aodEvent); - fHistPt->Fill(aodO->E()); - } - - - // if(esdEvent) { - // for(int icl = 0; icl < fESD->GetNumberOfCaloClusters(); icl++) { - // AliESDCaloCluster * cluster = fESD->GetCaloCluster(icl); - // if(cluster->GetEmcCpvDistance() > 10) { - // if (cluster->IsPHOS()) fHistPtPhos->Fill(cluster->E()); - // if (cluster->IsEMCAL()) fHistPtEmcal->Fill(cluster->E()); - // } - // } - // } - - PostData(1, fOutputList); } -//________________________________________________________________________ -void AliAnalysisTaskGammaJet::Terminate(Option_t *) -{ +//_____________________________________________________________________ +void AliAnalysisTaskGammaJet::Terminate(Option_t *) { // Draw result to the screen // Called once at the end of the query } -//________________________________________________________________________ +//_____________________________________________________________________ AliAODEvent * AliAnalysisTaskGammaJet::GetAODEvent() { //Get the AOD event from whereever it might be AliAODEvent * aodEvent = dynamic_cast(InputEvent()); if(!aodEvent) { aodEvent = AODEvent(); } - + return aodEvent; } -//________________________________________________________________________ -TClonesArray * AliAnalysisTaskGammaJet::GetConversionGammas() { +//_____________________________________________________________________ +TClonesArray * AliAnalysisTaskGammaJet::GetConversionGammas(const AliAODEvent * aodEvent) { - if( !(fDeltaAODFileName.Length() > 0) ) return NULL; + //Get Conversion gamma branch of AOD. First try standard AOD + TClonesArray * convGamma = dynamic_cast(aodEvent->FindListObject("GammaConv_gamma")); + + //If it's there, send it back + if(convGamma) return convGamma; + + //If AOD not in standard file have to locate it in delta AOD + if( !(fDeltaAODFileName.Length() > 0) ) return NULL; - //If AOD not in standard file have to locate it in extension AliAODHandler * aodHandler = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); if(aodHandler) { AliAODExtension * gExt = dynamic_cast(aodHandler->GetExtensions()->FindObject(fDeltaAODFileName)); if(gExt) { - AliAODEvent * aodEvent = gExt->GetAOD(); - return dynamic_cast(aodEvent->FindListObject("GammaConv")); + AliAODEvent * gcEvent = gExt->GetAOD(); + return dynamic_cast(gcEvent->FindListObject("GammaConv_gamma")); } } return NULL; } + + +//_____________________________________________________________________ +void AliAnalysisTaskGammaJet::FillPWG4PartCorrBranch( TClonesArray * gcBranch, TClonesArray * partCorrBranch , TString detector ) { + + for(int i = 0; i < gcBranch->GetEntriesFast(); i++) { + AliGammaConversionAODObject * gcObject = dynamic_cast(gcBranch->At(i)); + if ( gcObject ) { + AliAODPWG4ParticleCorrelation pc(gcObject->Px(), gcObject->Py(), gcObject->Pz(), gcObject->E()); + pc.SetTagged(gcObject->IsTagged()); + pc.SetTrackLabel(gcObject->GetLabel1(), gcObject->GetLabel2()); + pc.SetDetector(detector); + new((*partCorrBranch)[i]) AliAODPWG4ParticleCorrelation(pc); + + } else { + AliError(Form("Couldn't get gamma conversion aod object")); + } + + } +} + + +//_____________________________________________________________________ +AliAODPWG4ParticleCorrelation * AliAnalysisTaskGammaJet::PWG4PartFromGammaConvAODObject(AliGammaConversionAODObject * gcObject, TString detector ) { + + AliAODPWG4ParticleCorrelation * pc = new AliAODPWG4ParticleCorrelation(gcObject->Px(), gcObject->Py(), gcObject->Pz(), gcObject->E()); + pc->SetTagged(gcObject->IsTagged()); + pc->SetTrackLabel(gcObject->GetLabel1(), gcObject->GetLabel2()); + pc->SetDetector(detector); + return pc; +} + + +//_________________________________________________________________________ +void AliAnalysisTaskGammaJet::CleanUp() { + fPhotons->Delete(); +} + +//_________________________________________________________________________ +Bool_t AliAnalysisTaskGammaJet::IsIsolated( AliAODPWG4ParticleCorrelation * particle, TClonesArray * tracks, Float_t coneSize, Float_t ptThreshold ) { + //See header file for documentation + for(int it = 0; it < tracks->GetEntriesFast(); it++) { + if ( (it == particle->GetTrackLabel(0)) || it == particle->GetTrackLabel(1) ) + continue; + + //BALLE Svein:How are you checking the calorimeters for whether they are decay particles ? + + AliAODTrack * track = dynamic_cast(tracks->At(it)); + if (track) { + if ( IsInCone(particle->Eta() - track->Eta(), particle->Phi() - track->Phi(), coneSize) ) { + if (track->Pt() > ptThreshold) { + return kFALSE; + } + } + } else { + AliError(Form("Bad track!!!! ")); + } + } + + //No particle above threshold, it's isolated + return kTRUE; +} + + +//______________________________________________________________________________________________ +void AliAnalysisTaskGammaJet::ProcessCalorimeters( const AliAODEvent * const aodEvent ) { + + TClonesArray * clusters = aodEvent->GetCaloClusters(); + + + for(int ic = 0; ic < clusters->GetEntriesFast(); ic++) { + AliAODCaloCluster * cluster = dynamic_cast(clusters->At(ic)); + if (!cluster) { + AliError(Form("Error getting cluster")); + continue; + } + + + if (cluster->GetNCells() < 6) continue; + if (cluster->GetEmcCpvDistance() < 15) continue; + + TLorentzVector tlvec; + + AliAODVertex * vertex = aodEvent->GetPrimaryVertex(); + Double_t vertexPosition[3]; + vertex->GetXYZ(vertexPosition); + cluster->GetMomentum(tlvec, vertexPosition); + if (tlvec.Pt() < GetMinPt()) continue; + + AliAODPWG4ParticleCorrelation * photon = new AliAODPWG4ParticleCorrelation(tlvec); + + photon->SetIsolated( IsIsolated(photon, aodEvent->GetTracks(), GetConeSize(), GetPtThreshold()) ); + CorrelateWithJets(photon, aodEvent->GetJets()); + } + +} +//___________________________________________________________________________________________ +void AliAnalysisTaskGammaJet::ProcessConvGamma( const AliAODEvent * const aodEvent ) { + + TClonesArray * convGamma = GetConversionGammas(aodEvent); + if(!convGamma) { + AliError(Form("No convgamma")); + return; + } + + for (Int_t iPhot = 0; iPhot < convGamma->GetEntriesFast(); iPhot++) { + AliGammaConversionAODObject * aodO = dynamic_cast(convGamma->At(iPhot)); + + if (!aodO) { + AliError(Form("ERROR: Could not receive ga %d\n", iPhot)); + continue; + } + + if(aodO->Pt() < GetMinPt()) continue; + + + //Use the AODPWG4PartCorr shit! + AliAODPWG4ParticleCorrelation * photon = PWG4PartFromGammaConvAODObject(aodO, "ConvGamma"); + photon->SetIsolated( IsIsolated(photon, aodEvent->GetTracks(), GetConeSize(), GetPtThreshold()) ); + + + // if ( (aodO->Phi()) < 0 ) + // cout << aodO->Phi() << endl; + + CorrelateWithJets(photon, aodEvent->GetJets()); + + fHistPt->Fill(photon->Pt()); + delete photon; + + } +} + +void AliAnalysisTaskGammaJet::CorrelateWithJets(AliAODPWG4ParticleCorrelation * photon, const TClonesArray * const jets) { + //See header file for documentation + if (jets) { + for(int ij = 0; ij < jets->GetEntriesFast(); ij++) { + AliAODJet * jet = dynamic_cast(jets->At(ij)); + if(jet) { + fHistPtJets->Fill(jet->Pt()); + + Float_t dPhi = TMath::Abs(photon->Phi() - jet->Phi()); + if (photon->IsIsolated()) + fHistGammaJetsIso->Fill(dPhi, jet->Pt()/photon->Pt()); + else + fHistGammaJets->Fill(dPhi); + + } + } + } +} diff --git a/PWG4/GammaConv/AliAnalysisTaskGammaJet.h b/PWG4/GammaConv/AliAnalysisTaskGammaJet.h index 4fb26809ba7..936ec9079c0 100644 --- a/PWG4/GammaConv/AliAnalysisTaskGammaJet.h +++ b/PWG4/GammaConv/AliAnalysisTaskGammaJet.h @@ -6,6 +6,10 @@ class TH1F; class AliESDEvent; +class AliGammaConversionAODObject; +class AliAODPWG4ParticleCorrelation; +class TClonesArray; +class TString; #include "AliAnalysisTaskSE.h" @@ -20,23 +24,69 @@ class AliAnalysisTaskGammaJet : public AliAnalysisTaskSE { virtual void Terminate(Option_t *); void SetDeltaAODFileName(TString string) { fDeltaAODFileName = string;} + //Move all of these somewhere + inline Float_t GetMinPt() const { return fMinPt;} + void SetMinPt( const Float_t pt ) { fMinPt = pt; } + inline Float_t GetConeSize () const { return fConeSize; } + void SetConeSize ( const Float_t cs ) { fConeSize = cs; } + inline Float_t GetPtThreshold () const { return fPtThreshold; } + void SetPtThreshold ( Float_t ptt ) { fPtThreshold = ptt; } + + + + private: + //Clean up + void CleanUp(); + //Get the AOD event from whereever it might be accessible AliAODEvent * GetAODEvent(); - TClonesArray * GetConversionGammas(); + //Get Conversion gammas branch + TClonesArray * GetConversionGammas(const AliAODEvent * aodEvent); + + //Create PWG4 particles from the aod objects + AliAODPWG4ParticleCorrelation * PWG4PartFromGammaConvAODObject(AliGammaConversionAODObject * gcObject, TString detector); + + //Fill AOD tree with PWG4 particles + void FillPWG4PartCorrBranch( TClonesArray * gcBranch, TClonesArray * partCorrBranch, TString detector); + //Is particle isolated + Bool_t IsIsolated( AliAODPWG4ParticleCorrelation * particle, TClonesArray * tracks, Float_t coneSize, Float_t ptThreshold); + + //Process conv gamma + void ProcessConvGamma( const AliAODEvent * const aodEvent ); + + //Process calorimeters + void ProcessCalorimeters( const AliAODEvent * const aodEvent ); + + //Correlate particle with jets + void CorrelateWithJets(AliAODPWG4ParticleCorrelation * photon, const TClonesArray * const jets); + + //Is eta - phi distance smaller than conesize ? + inline Bool_t IsInCone(Float_t dEta, Float_t dPhi, Float_t coneSize) { + return ( (dEta*dEta + dPhi*dPhi) < coneSize*coneSize); + } + TList *fOutputList; //! Output list TH1F *fHistPt; //! Pt spectrum TH1F *fHistPtPhos; //! Pt spectrum TH1F *fHistPtEmcal; //! Pt spectrum TH1F *fHistPtJets; //! Pt spectrum TH1F *fHistGammaJets; //!Phi correlations - + TH1F *fHistGammaJetsIso; //!Phi correlations + + + Float_t fMinPt; //Minimum pt for correlation + Float_t fConeSize; //cone size for isolation + Float_t fPtThreshold; //Threshold pt for isolation + TString fDeltaAODFileName;//! File where Gamma Conv AOD is located, if not in default AOD + TClonesArray * fPhotons; + AliAnalysisTaskGammaJet(const AliAnalysisTaskGammaJet&); // not implemented AliAnalysisTaskGammaJet& operator=(const AliAnalysisTaskGammaJet&); // not implemented -- 2.43.5