From: kleinb Date: Tue, 28 Sep 2010 07:35:25 +0000 (+0000) Subject: Adding the background calculation to the cluster task (Leticia) X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=d6e66a820c40efca661d2c7b28ec1dfa4637b7ff Adding the background calculation to the cluster task (Leticia) --- diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index 36cf719c03c..161eced5f6a 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -55,7 +55,7 @@ #include "AliJetKineReaderHeader.h" #include "AliGenCocktailEventHeader.h" #include "AliInputEventHandler.h" - +#include "AliAODJetEventBackground.h" #include "fastjet/PseudoJet.hh" #include "fastjet/ClusterSequenceArea.hh" @@ -94,6 +94,9 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fStrategy(fastjet::Best), fRecombScheme(fastjet::BIpt_scheme), fAreaType(fastjet::active_area), + fGhostArea(0.01), + fActiveAreaRepeats(1), + fGhostEtamax(1.5), fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), @@ -171,6 +174,9 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name): fStrategy(fastjet::Best), fRecombScheme(fastjet::BIpt_scheme), fAreaType(fastjet::active_area), + fGhostArea(0.01), + fActiveAreaRepeats(1), + fGhostEtamax(1.5), fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), @@ -250,7 +256,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n"); - + if(fNonStdBranch.Length()!=0) { @@ -262,6 +268,20 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() TClonesArray *tca = new TClonesArray("AliAODJet", 0); tca->SetName(fNonStdBranch.Data()); AddAODBranch("TClonesArray",&tca,fNonStdFile.Data()); + + + TClonesArray *tcaran = new TClonesArray("AliAODJet", 0); + tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); + AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data()); + + + if(!AODEvent() || !AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){ + AliAODJetEventBackground* evBkg = new AliAODJetEventBackground(); + evBkg->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data())); + AddAODBranch("AliAODJetEventBackground",&evBkg,fNonStdFile.Data()); + + } + if(fNonStdFile.Length()!=0){ // // case that we have an AOD extension we need to fetch the jets from the extended output @@ -509,12 +529,24 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // handle and reset the output jet branch // only need this once TClonesArray* jarray = 0; + TClonesArray* jarrayran = 0; + AliAODJetEventBackground* evBkg = 0; if(fNonStdBranch.Length()!=0) { if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data())); if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data())); if(jarray)jarray->Delete(); // this is our responsibility, clear before filling again + if(AODEvent())jarrayran = (TClonesArray*)(AODEvent()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random"))); + if(!jarrayran)jarrayran = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",fNonStdBranch.Data(),"random"))); + if(jarrayran)jarrayran->Delete(); // this is our responsibility, clear before filling again + + if(AODEvent())evBkg = (AliAODJetEventBackground*)(AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))); + if(!evBkg) evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))); + if(evBkg)evBkg->Reset(); + } + + // // Execute analysis for current event // @@ -630,19 +662,25 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // run fast jet // employ setters for these... - double ghostEtamax = 0.9; - double ghostArea = 0.01; - int activeAreaRepeats = 1; + // now create the object that holds info about ghosts - fastjet::GhostedAreaSpec ghostSpec(ghostEtamax, activeAreaRepeats, ghostArea)\ -; + fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea); fastjet::AreaType areaType = fastjet::active_area; fastjet::AreaDefinition areaDef = fastjet::AreaDefinition(areaType,ghostSpec); - - fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy); + fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy); vector inclusiveJets, sortedJets; fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef); + //range where to compute background + Double_t phiMin = 0, phiMax = 0, rapMin = 0, rapMax = 0; + phiMin = 0; + phiMax = 2*TMath::Pi(); + rapMax = fGhostEtamax - fRparam; + rapMin = - fGhostEtamax + fRparam; + fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax); + + + inclusiveJets = clustSeq.inclusive_jets(); sortedJets = sorted_by_pt(inclusiveJets); @@ -759,6 +797,27 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); } delete recIter; + + //background estimates:all bckg jets(0) & wo the 2 hardest(1) + + if(evBkg){ + vector jets2=sortedJets; + if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); + Double_t bkg1=0; + Double_t sigma1=0.; + Double_t meanarea1=0.; + Double_t bkg2=0; + Double_t sigma2=0.; + Double_t meanarea2=0.; + + clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, false); + evBkg->SetBackground(0,bkg1,sigma1,meanarea1); + clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, false); + evBkg->SetBackground(1,bkg2,sigma2,meanarea2); + } + + + } // fill track information @@ -815,7 +874,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // find the random jets vector inclusiveJetsRan, sortedJetsRan; - fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef); + fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); inclusiveJetsRan = clustSeqRan.inclusive_jets(); sortedJetsRan = sorted_by_pt(inclusiveJetsRan); @@ -864,8 +923,8 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh2TracksLeadingJetPhiPtWRan->Fill(dPhi,pt,tmpPt); } - - + Int_t nAodOutJetsRan = 0; + AliAODJet *aodOutJetRan = 0; for(int j = 0; j < nRecRan;j++){ AliAODJet tmpRec (sortedJetsRan[j].px(), sortedJetsRan[j].py(), sortedJetsRan[j].pz(), sortedJetsRan[j].E()); Float_t tmpPt = tmpRec.Pt(); @@ -876,6 +935,27 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh2NConstPtRan->Fill(tmpPt,constituents.size()); fh2PtNchRan->Fill(nCh,tmpPt); fh2PtNchNRan->Fill(nCh,tmpPt,constituents.size()); + + + if(tmpPt>fJetOutputMinPt&&jarrayran){ + aodOutJetRan = new ((*jarrayran)[nAodOutJetsRan++]) AliAODJet(tmpRec); + Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); + + aodOutJetRan->SetEffArea(arearan,0); } + + + + + for(unsigned int ic = 0; ic < constituents.size();ic++){ + // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); + // fh1PtJetConstRec->Fill(part->Pt()); + if(aodOutJetRan){ + aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index())); + }} + + + + // correlation Float_t tmpPhi = tmpRec.Phi(); if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; @@ -887,6 +967,18 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) continue; } } + + + if(evBkg){ + Double_t bkg3=0.; + Double_t sigma3=0.; + Double_t meanarea3=0.; + clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, false); + evBkg->SetBackground(2,bkg3,sigma3,meanarea3); + } + + + } diff --git a/JETAN/AliAnalysisTaskJetCluster.h b/JETAN/AliAnalysisTaskJetCluster.h index 70cb809af9e..4c8722c5ebc 100644 --- a/JETAN/AliAnalysisTaskJetCluster.h +++ b/JETAN/AliAnalysisTaskJetCluster.h @@ -23,7 +23,8 @@ class AliAODExtension; class AliAODJet; class AliGenPythiaEventHeader; class AliCFManager; - +class AliAODJetEventBackground; +class AliJetFinder; class TList; class TChain; class TH2F; @@ -47,6 +48,8 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual void Terminate(Option_t *option); virtual Bool_t Notify(); + + virtual void SetUseGlobalSelection(Bool_t b){fUseGlobalSelection = b;} virtual void SetAODTrackInput(Bool_t b){fUseAODTrackInput = b;} virtual void SetAODMCInput(Bool_t b){fUseAODMCInput = b;} @@ -55,7 +58,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE virtual void SetTrackTypeRec(Int_t i){fTrackTypeRec = i;} virtual void SetTrackPtCut(Float_t x){fTrackPtCut = x;} virtual void SetFilterMask(UInt_t i){fFilterMask = i;} - + virtual void SetNSkipLeadingRan(Int_t x){fNSkipLeadingRan = x;} virtual void SetJetOutputBranch(const char *c){fNonStdBranch = c;} @@ -73,6 +76,11 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE void SetStrategy(fastjet::Strategy f) {fStrategy = f;} void SetRecombScheme(fastjet::RecombinationScheme f) {fRecombScheme = f;} void SetAreaType(fastjet::AreaType f) {fAreaType = f;} + void SetGhostArea(Double_t f) {fGhostArea = f;} + void SetActiveAreaRepeats(Int_t f) {fActiveAreaRepeats = f;} + void SetGhostEtamax(Double_t f) {fGhostEtamax = f;} + + // Helper // @@ -121,64 +129,66 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE fastjet::Strategy fStrategy; //= fastjet::Best; fastjet::RecombinationScheme fRecombScheme; // = fastjet::BIpt_scheme; fastjet::AreaType fAreaType; - - TProfile* fh1Xsec; // pythia cross section and trials - TH1F* fh1Trials; // trials are added - TH1F* fh1PtHard; // Pt har of the event... - TH1F* fh1PtHardNoW; // Pt har of the event without weigt - TH1F* fh1PtHardTrials; // Number of trials - - TH1F* fh1NJetsRec; // number of reconstructed jets - TH1F* fh1NConstRec;// number of constiutens in leading jet - TH1F* fh1NConstLeadingRec;// number of constiutens in leading jet - TH1F* fh1PtJetsRecIn; // Jet pt for all jets - TH1F* fh1PtJetsLeadingRecIn; // Jet pt for all jets - TH1F* fh1PtJetConstRec;// pt of constituents + Double_t fGhostArea; + Int_t fActiveAreaRepeats; + Double_t fGhostEtamax; + TProfile* fh1Xsec; //! pythia cross section and trials + TH1F* fh1Trials; //! trials are added + TH1F* fh1PtHard; //! Pt har of the event... + TH1F* fh1PtHardNoW; //! Pt har of the event without weigt + TH1F* fh1PtHardTrials; //! Number of trials + + TH1F* fh1NJetsRec; //! number of reconstructed jets + TH1F* fh1NConstRec;//! number of constiutens in leading jet + TH1F* fh1NConstLeadingRec;//! number of constiutens in leading jet + TH1F* fh1PtJetsRecIn; //! Jet pt for all jets + TH1F* fh1PtJetsLeadingRecIn; //! Jet pt for all jets + TH1F* fh1PtJetConstRec;//! pt of constituents TH1F* fh1PtJetConstLeadingRec;// pt of constituents - TH1F* fh1PtTracksRecIn; // track pt for all tracks - TH1F* fh1PtTracksLeadingRecIn; // track pt for all tracks + TH1F* fh1PtTracksRecIn; //! track pt for all tracks + TH1F* fh1PtTracksLeadingRecIn; //! track pt for all tracks // Randomized track histos - TH1F* fh1NJetsRecRan; // number of reconstructed jets from randomized - TH1F* fh1NConstRecRan;// number of constiutens in leading jet - TH1F* fh1PtJetsLeadingRecInRan; // Jet pt for all jets - TH1F* fh1NConstLeadingRecRan;// number of constiutens in leading jet - TH1F* fh1PtJetsRecInRan; // Jet pt for all jets - - TH1F* fh1PtTracksGenIn; // track pt for all tracks - TH1F* fh1Nch; // charged particle mult - - TH2F* fh2NRecJetsPt; // Number of found jets above threshold - TH2F* fh2NRecTracksPt; // Number of found tracks above threshold - TH2F* fh2NConstPt; // number of constituents vs. pt - TH2F* fh2NConstLeadingPt; // number of constituents vs. pt - TH2F* fh2JetPhiEta; // jet phi eta - TH2F* fh2LeadingJetPhiEta; // leading jet phi eta - TH2F* fh2JetEtaPt; // leading jet eta - TH2F* fh2LeadingJetEtaPt; // leading jet eta - TH2F* fh2TrackEtaPt; // track eta - TH2F* fh2LeadingTrackEtaPt; // leading track eta - TH2F* fh2JetsLeadingPhiEta; // jet phi eta - TH2F* fh2JetsLeadingPhiPt; // jet correlation with leading jet - TH2F* fh2TracksLeadingPhiEta; // track correlation with leading track - TH2F* fh2TracksLeadingPhiPt; // track correlation with leading track - TH2F* fh2TracksLeadingJetPhiPt; // track correlation with leading Jet - TH2F* fh2JetsLeadingPhiPtW; // jet correlation with leading jet - TH2F* fh2TracksLeadingPhiPtW; // track correlation with leading track - TH2F* fh2TracksLeadingJetPhiPtW; // track correlation with leading Jet - TH2F* fh2NRecJetsPtRan; // Number of found jets above threshold - TH2F* fh2NConstPtRan; // number of constituents vs. pt - TH2F* fh2NConstLeadingPtRan; // number of constituents vs. pt - TH2F* fh2PtNch; // p_T of cluster vs. multiplicity, - TH2F* fh2PtNchRan; // p_T of cluster vs. multiplicity,random - TH2F* fh2PtNchN; // p_T of cluster vs. multiplicity, weigthed with constituents - TH2F* fh2PtNchNRan; // p_T of cluster vs. multiplicity, weigthed with constituents random - TH2F* fh2TracksLeadingJetPhiPtRan; // track correlation with leading Jet - TH2F* fh2TracksLeadingJetPhiPtWRan; // track correlation with leading Jet + TH1F* fh1NJetsRecRan; //! number of reconstructed jets from randomized + TH1F* fh1NConstRecRan;//! number of constiutens in leading jet + TH1F* fh1PtJetsLeadingRecInRan; //! Jet pt for all jets + TH1F* fh1NConstLeadingRecRan;//! number of constiutens in leading jet + TH1F* fh1PtJetsRecInRan; //! Jet pt for all jets + + TH1F* fh1PtTracksGenIn; //! track pt for all tracks + TH1F* fh1Nch; //! charged particle mult + + TH2F* fh2NRecJetsPt; //! Number of found jets above threshold + TH2F* fh2NRecTracksPt; //! Number of found tracks above threshold + TH2F* fh2NConstPt; //! number of constituents vs. pt + TH2F* fh2NConstLeadingPt; //! number of constituents vs. pt + TH2F* fh2JetPhiEta; //! jet phi eta + TH2F* fh2LeadingJetPhiEta; //! leading jet phi eta + TH2F* fh2JetEtaPt; //! leading jet eta + TH2F* fh2LeadingJetEtaPt; //! leading jet eta + TH2F* fh2TrackEtaPt; //! track eta + TH2F* fh2LeadingTrackEtaPt; //! leading track eta + TH2F* fh2JetsLeadingPhiEta; //! jet phi eta + TH2F* fh2JetsLeadingPhiPt; //! jet correlation with leading jet + TH2F* fh2TracksLeadingPhiEta; //! track correlation with leading track + TH2F* fh2TracksLeadingPhiPt; //! track correlation with leading track + TH2F* fh2TracksLeadingJetPhiPt; //! track correlation with leading Jet + TH2F* fh2JetsLeadingPhiPtW; //! jet correlation with leading jet + TH2F* fh2TracksLeadingPhiPtW; //! track correlation with leading track + TH2F* fh2TracksLeadingJetPhiPtW; //! track correlation with leading Jet + TH2F* fh2NRecJetsPtRan; //! Number of found jets above threshold + TH2F* fh2NConstPtRan; //! number of constituents vs. pt + TH2F* fh2NConstLeadingPtRan; //! number of constituents vs. pt + TH2F* fh2PtNch; //! p_T of cluster vs. multiplicity, + TH2F* fh2PtNchRan; //! p_T of cluster vs. multiplicity,random + TH2F* fh2PtNchN; //! p_T of cluster vs. multiplicity, weigthed with constituents + TH2F* fh2PtNchNRan; //! p_T of cluster vs. multiplicity, weigthed with constituents random + TH2F* fh2TracksLeadingJetPhiPtRan; //! track correlation with leading Jet + TH2F* fh2TracksLeadingJetPhiPtWRan; //! track correlation with leading Jet TList *fHistList; // Output list - ClassDef(AliAnalysisTaskJetCluster, 4) + ClassDef(AliAnalysisTaskJetCluster, 5) }; #endif diff --git a/PWG4/macros/AddTaskJetCluster.C b/PWG4/macros/AddTaskJetCluster.C index 2638bbffaba..f65beb5a346 100644 --- a/PWG4/macros/AddTaskJetCluster.C +++ b/PWG4/macros/AddTaskJetCluster.C @@ -123,7 +123,7 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte if(TMath::Abs(radius-0.4)<0.01){ pwg4spec->SetJetOutputBranch(Form("clusters%s_%s%s",bRec,jf,cRadius)); - pwg4spec->SetJetOutputMinPt(1); // store only jets / clusters above a certain threshold + pwg4spec->SetJetOutputMinPt(0); // store only jets / clusters above a certain threshold } pwg4spec->SetNSkipLeadingRan(nSkip);