#include "AliAnalysisTaskFastEmbedding.h"
#include "AliAnalysisManager.h"
+#include "AliESDtrack.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
,fEvtSelecMode(0)
,fEvtSelMinJetPt(-1)
,fEvtSelMaxJetPt(-1)
+ ,fEvtSelMinJetEta(-999.)
+ ,fEvtSelMaxJetEta( 999.)
+ ,fEvtSelMinJetPhi(0.)
+ ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
,fToyMinNbOfTracks(1)
,fToyMaxNbOfTracks(1)
,fToyMinTrackPt(50.)
,fh1TrackPt(0)
,fh2TrackEtaPhi(0)
,fh1TrackN(0)
+ ,fh1JetPt(0)
+ ,fh2JetEtaPhi(0)
+ ,fh1JetN(0)
,fh1MCTrackPt(0)
,fh2MCTrackEtaPhi(0)
,fh1MCTrackN(0)
+ ,fh1AODfile(0)
{
,fEvtSelecMode(0)
,fEvtSelMinJetPt(-1)
,fEvtSelMaxJetPt(-1)
+ ,fEvtSelMinJetEta(-999.)
+ ,fEvtSelMaxJetEta( 999.)
+ ,fEvtSelMinJetPhi(0.)
+ ,fEvtSelMaxJetPhi(TMath::Pi()*2.)
,fToyMinNbOfTracks(1)
,fToyMaxNbOfTracks(1)
,fToyMinTrackPt(50.)
,fh1TrackPt(0)
,fh2TrackEtaPhi(0)
,fh1TrackN(0)
+ ,fh1JetPt(0)
+ ,fh2JetEtaPhi(0)
+ ,fh1JetN(0)
,fh1MCTrackPt(0)
,fh2MCTrackEtaPhi(0)
,fh1MCTrackN(0)
+ ,fh1AODfile(0)
{
// constructor
DefineOutput(1, TList::Class());
,fEvtSelecMode(copy.fEvtSelecMode)
,fEvtSelMinJetPt(copy.fEvtSelMinJetPt)
,fEvtSelMaxJetPt(copy.fEvtSelMaxJetPt)
+ ,fEvtSelMinJetEta(copy.fEvtSelMinJetEta)
+ ,fEvtSelMaxJetEta(copy.fEvtSelMaxJetEta)
+ ,fEvtSelMinJetPhi(copy.fEvtSelMinJetPhi)
+ ,fEvtSelMaxJetPhi(copy.fEvtSelMaxJetPhi)
,fToyMinNbOfTracks(copy.fToyMinNbOfTracks)
,fToyMaxNbOfTracks(copy.fToyMaxNbOfTracks)
,fToyMinTrackPt(copy.fToyMinTrackPt)
,fh1TrackPt(copy.fh1TrackPt)
,fh2TrackEtaPhi(copy.fh2TrackEtaPhi)
,fh1TrackN(copy.fh1TrackN)
+ ,fh1JetPt(copy.fh1JetPt)
+ ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
+ ,fh1JetN(copy.fh1JetN)
,fh1MCTrackPt(copy.fh1MCTrackPt)
,fh2MCTrackEtaPhi(copy.fh2MCTrackEtaPhi)
,fh1MCTrackN(copy.fh1MCTrackN)
+ ,fh1AODfile(copy.fh1AODfile)
{
// copy constructor
}
fEvtSelecMode = o.fEvtSelecMode;
fEvtSelMinJetPt = o.fEvtSelMinJetPt;
fEvtSelMaxJetPt = o.fEvtSelMaxJetPt;
+ fEvtSelMinJetEta = o.fEvtSelMinJetEta;
+ fEvtSelMaxJetEta = o.fEvtSelMaxJetEta;
+ fEvtSelMinJetPhi = o.fEvtSelMinJetPhi;
+ fEvtSelMaxJetPhi = o.fEvtSelMaxJetPhi;
fToyMinNbOfTracks = o.fToyMinNbOfTracks;
fToyMaxNbOfTracks = o.fToyMaxNbOfTracks;
fToyMinTrackPt = o.fToyMinTrackPt;
fh1TrackPt = o.fh1TrackPt;
fh2TrackEtaPhi = o.fh2TrackEtaPhi;
fh1TrackN = o.fh1TrackN;
+ fh1JetPt = o.fh1JetPt;
+ fh2JetEtaPhi = o.fh2JetEtaPhi;
+ fh1JetN = o.fh1JetN;
fh1MCTrackPt = o.fh1MCTrackPt;
fh2MCTrackEtaPhi = o.fh2MCTrackEtaPhi;
fh1MCTrackN = o.fh1MCTrackN;
+ fh1AODfile = o.fh1AODfile;
}
return *this;
// create output objects
if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::UserCreateOutputObjects()");
- rndm = new TRandom3();
- Int_t id = GetJobID();
- if(id>-1) rndm->SetSeed(id);
- else rndm->SetSeed(); // a TTUID is generated and used for seed
- AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
-
-
- // embed mode with AOD
- if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-
- // open input AOD
- if(fAODPathArray){
- Int_t rc = SelectAODfile();
- if(rc<0) return;
- }
-
- Int_t rc = OpenAODfile();
- if(rc<0) return;
-
- } //end: embed mode with AOD
-
-
// connect output aod
// create a new branch for extra tracks
fAODout = AODEvent();
//qa histograms
OpenFile(1);
- fHistList = new TList();
+ if(!fHistList) fHistList = new TList();
+ fHistList->SetOwner(kTRUE);
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
- fh1TrackPt = new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 300, 0., 300.);
+ fh1TrackPt = new TH1F("fh1TrackPt","pT of extra tracks;p_{T};entries", 250, 0., 250.);
fh2TrackEtaPhi = new TH2F("fh2TrackEtaPhi","eta-phi distribution of extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
- fh1TrackN = new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 3000.);
+ fh1TrackN = new TH1F("fh1TrackN", "nb. of extra tracks per event;nb. of tracks;entries",300, 0., 300.);
fHistList->Add(fh1TrackPt);
fHistList->Add(fh2TrackEtaPhi);
fHistList->Add(fh1TrackN);
+
+ if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+
+ fh1JetPt = new TH1F("fh1JetPt", "pT of extra jets;p_{T};entries", 250, 0., 250.);
+ fh2JetEtaPhi = new TH2F("fh2JetEtaPhi", "eta-phi distribution of extra jets;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
+ fh1JetN = new TH1F("fh1JetN", "nb. of extra jets per event;nb. of jets;entries",20,0.,20.);
+
+ fHistList->Add(fh1JetPt);
+ fHistList->Add(fh2JetEtaPhi);
+ fHistList->Add(fh1JetN);
+ }
if(fAODevent && fAODevent->FindListObject("mcparticles") && strlen(fMCparticlesBranch.Data())){
- fh1MCTrackPt = new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 300, 0., 300.);
+ fh1MCTrackPt = new TH1F("fh1MCTrackPt","pT of MC extra tracks;p_{T};entries", 250, 0., 250.);
fh2MCTrackEtaPhi = new TH2F("fh2MCTrackEtaPhi","eta-phi distribution of MC extra tracks;#eta;#phi", 20, -1., 1., 60, 0., 2*TMath::Pi());
- fh1MCTrackN = new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 3000.);
-
+ fh1MCTrackN = new TH1F("fh1MCTrackN", "nb. of MC extra tracks per event;nb. of tracks;entries",300, 0., 300.);
+
fHistList->Add(fh1MCTrackPt);
fHistList->Add(fh2MCTrackEtaPhi);
fHistList->Add(fh1MCTrackN);
-
+
+ }
+
+ fh1AODfile = new TH1I("fh1AODfile", "overview of opened AOD files from the array", 100, 0, 100);
+ fHistList->Add(fh1AODfile);
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
+ TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
+ if (h1){
+ h1->Sumw2();
+ continue;
+ }
}
TH1::AddDirectory(oldStatus);
+
+ // set seed
+ rndm = new TRandom3();
+ Int_t id = GetJobID();
+ if(id>-1) rndm->SetSeed(id);
+ else rndm->SetSeed(); // a TTUID is generated and used for seed
+ AliInfo(Form("TRandom3 seed: %d", rndm->GetSeed()));
+
+ // embed mode with AOD
+ if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
+
+ // open input AOD
+ Int_t rc = OpenAODfile();
+ if(rc<0) return;
+ fh1AODfile->Fill(rc);
+
+ } //end: embed mode with AOD
+
+
+
+ PostData(1, fHistList);
}
if(!fAODout){
AliError("Need output AOD, but is not connected.");
+ PostData(1, fHistList);
return;
}
TClonesArray *tracks = (TClonesArray*)(fAODout->FindListObject(fTrackBranch.Data()));
if(!tracks){
AliError("Extra track branch not found in output.");
+ PostData(1, fHistList);
return;
}
tracks->Delete();
if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
if(!fAODevent){
AliError("Need input AOD, but is not connected.");
+ PostData(1, fHistList);
return;
}
else aodJets = fAODevent->GetJets();
if(!aodJets){
AliError("Could not find jets in AOD. Check jet branch when indicated.");
+ PostData(1, fHistList);
return;
}
}
else {
AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
- Int_t rc = SelectAODfile();
- if(rc<0) return;
- rc = OpenAODfile();
- if(rc<0) return;
+ Int_t rc = OpenAODfile();
+ if(rc<0) {
+ PostData(1, fHistList);
+ return;
+ }
+ fh1AODfile->Fill(rc);
// new file => we must use the new jet array
if(fJetBranch.Length()) aodJets = dynamic_cast<TClonesArray*>(fAODevent->FindListObject(fJetBranch.Data()));
else aodJets = fAODevent->GetJets();
if(!aodJets){
AliError("Could not find jets in AOD. Check jet branch when indicated.");
+ PostData(1, fHistList);
return;
}
}
if(!jet) continue;
if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
- && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)){
+ && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+ && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+ && (jet->Phi()>=fEvtSelMinJetPhi && jet->Phi()<=fEvtSelMaxJetPhi)){
useEntry = kTRUE;
break;
}
if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks){ // take all tracks or jet tracks
// loop over input tracks
// add to output aod
- Int_t nTracks = fAODevent->GetNumberOfTracks();
+ Int_t nTracks = 0;
+ Int_t nJets = aodJets->GetEntries();
+ Int_t nSelectedJets = 0;
+ AliAODJet *leadJet = 0x0; // used for jet tracks only
+
+ if(fEmbedMode==kAODFull){
+ nTracks = fAODevent->GetNumberOfTracks();
+
+ for(Int_t ij=0; ij<nJets; ++ij){
+ AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+ if(!jet) continue;
+ if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+ && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+ && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+ && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
+
+ fh1JetPt->Fill(jet->Pt());
+ fh2JetEtaPhi->Fill(jet->Eta(), jet->Phi());
+ nSelectedJets++;
+
+ }
+ }
+ fh1JetN->Fill(nSelectedJets);
+ }
+
+ if(fEmbedMode==kAODJetTracks){
+ // look for leading jet within selection
+ // get reference tracks for this jet
+ Float_t leadJetPt = 0.;
+ for(Int_t ij=0; ij<nJets; ++ij){
+ AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+ if(!jet) continue;
+ if( (fEvtSelMinJetPt==-1. || jet->Pt()>=fEvtSelMinJetPt)
+ && (fEvtSelMaxJetPt==-1. || jet->Pt()<=fEvtSelMaxJetPt)
+ && (jet->Eta()>=fEvtSelMinJetEta && jet->Eta()<=fEvtSelMaxJetEta)
+ && (jet->Phi()>=fEvtSelMinJetPhi && jet->Eta()<=fEvtSelMaxJetPhi)){
+ if(jet->Pt()>leadJetPt){
+ leadJetPt = jet->Pt();
+ leadJet = jet;
+ }
+ }
+ }
+ if(leadJet){
+ nTracks = leadJet->GetRefTracks()->GetEntriesFast();
+ fh1JetPt->Fill(leadJet->Pt());
+ fh2JetEtaPhi->Fill(leadJet->Eta(), leadJet->Pt());
+ fh1JetN->Fill(1);
+ }
+ }
+
fh1TrackN->Fill((Float_t)nTracks);
for(Int_t it=0; it<nTracks; ++it){
- AliAODTrack *tmpTr = fAODevent->GetTrack(it);
- // ?? test filter bit, or something else??
+ AliAODTrack *tmpTr = 0x0;
+ if(fEmbedMode==kAODFull) tmpTr = fAODevent->GetTrack(it);
+ if(fEmbedMode==kAODJetTracks) tmpTr = dynamic_cast<AliAODTrack*>(leadJet->GetRefTracks()->At(it));
+ if(!tmpTr) continue;
- if(fEmbedMode==kAODJetTracks){
- // jet track? else continue
- // not implemented yet
- continue;
- }
+ tmpTr->SetStatus(AliESDtrack::kEmbedded);
new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
dummy = (*tracks)[nAODtracks-1];
AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
if(!jet) continue;
AliAODTrack *tmpTr = (AliAODTrack*)jet;
+ tmpTr->SetFlags(AliESDtrack::kEmbedded);
new ((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
dummy = (*tracks)[nAODtracks-1];
Double_t phi = rndm->Uniform(fToyMinTrackPhi,fToyMaxTrackPhi);
phi = TVector2::Phi_0_2pi(phi);
- Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
+ if(fDebug) Printf("Add track #%d: pt %.2f, eta %.2f, phi %.2f", i, pt, eta, phi);
Double_t theta = 2*TMath::ATan(1/TMath::Exp(eta));
Float_t mom[3];
fToyFilterMap, // select info
-999. // chi2 per NDF
);
- tmpTr->SetFlags(1<<27);
+ tmpTr->SetFlags(AliESDtrack::kEmbedded);
new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
dummy = (*tracks)[nAODtracks-1];
return n;
}
-Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(){
+//__________________________________________________________________________
+
+Int_t AliAnalysisTaskFastEmbedding::OpenAODfile(Int_t trial){
+
+ if(trial>3){
+ AliError("Could not open AOD files, give up ...");
+ return -1;
+ }
+
+ Int_t rc = 0;
+ if(fAODPathArray) rc = SelectAODfile();
+ if(rc<0) return -1;
TDirectory *owd = gDirectory;
if (fAODfile)
fAODfile = TFile::Open(fAODPath.Data());
owd->cd();
if(!fAODfile){
- AliError("Could not open AOD file.");
- return -1;
+
+ rc = -1;
+ if(fAODPathArray){
+ AliError(Form("Could not open AOD file %s, try another from the list ...", fAODPath.Data()));
+ rc = OpenAODfile(trial+1);
+ } else {
+ AliError(Form("Could not open AOD file %s.", fAODPath.Data()));
+ }
+
+ return rc;
}
fAODtree = (TTree*)fAODfile->Get("aodTree");
delete fAODevent;
fAODevent = new AliAODEvent();
fAODevent->ReadFromTree(fAODtree);
- return 1;
+ return rc; // file position in AOD path array, if array available
}
#include "TMath.h"
#include "TH1F.h"
#include "TH2F.h"
+#include "TH3F.h"
#include "TCanvas.h"
#include "AliLog.h"
fMinContribVtx(1),
fVtxZMin(-8.),
fVtxZMax(8.),
- fEvtClassMin(0),
- fEvtClassMax(10),
+ fEvtClassMin(1),
+ fEvtClassMax(4),
fCentMin(0.),
fCentMax(100.),
fJetEtaMin(-.5),
fJetEtaMax(.5),
- fJetDeltaEta(.4),
- fJetDeltaPhi(.4),
+ fJetPtMin(20.),
+ fJetPtFractionMin(0.5),
+ fNMatchJets(4),
+ //fJetDeltaEta(.4),
+ //fJetDeltaPhi(.4),
fkNbranches(2),
- fkEvtClasses(10),
+ fkEvtClasses(5),
fOutputList(0x0),
fHistEvtSelection(0x0),
- fHistPtLeadingJet(new TH1F*[fkNbranches]),
- fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
- fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
- fHistDeltaEtaDeltaPhiLeadingJet(0x0),
- fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
- fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
+ fHistEvtClass(0x0),
+ fHistCentrality(0x0),
+ fHistPtJet(new TH1F*[fkNbranches]),
+ fHistEtaPhiJet(new TH2F*[fkNbranches]),
+ fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
+ fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
+ fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
+ fHistPtFraction(new TH2F*[fkEvtClasses]),
fHistPtPtExtra(0x0),
fHistPtResponse(new TH2F*[fkEvtClasses]),
fHistPtSmearing(new TH2F*[fkEvtClasses]),
- fHistdR(new TH2F*[fkEvtClasses]),
- fHistArea(new TH2F*[fkEvtClasses])
+ fHistDeltaR(new TH2F*[fkEvtClasses]),
+ fHistArea(new TH2F*[fkEvtClasses]),
+ fHistDeltaArea(new TH2F*[fkEvtClasses])
{
// default Constructor
fMinContribVtx(1),
fVtxZMin(-8.),
fVtxZMax(8.),
- fEvtClassMin(0),
- fEvtClassMax(10),
+ fEvtClassMin(1),
+ fEvtClassMax(4),
fCentMin(0.),
fCentMax(100.),
fJetEtaMin(-.5),
fJetEtaMax(.5),
- fJetDeltaEta(.4),
- fJetDeltaPhi(.4),
+ fJetPtMin(20.),
+ fJetPtFractionMin(0.5),
+ fNMatchJets(4),
+ //fJetDeltaEta(.4),
+ //fJetDeltaPhi(.4),
fkNbranches(2),
- fkEvtClasses(10),
+ fkEvtClasses(5),
fOutputList(0x0),
fHistEvtSelection(0x0),
- fHistPtLeadingJet(new TH1F*[fkNbranches]),
- fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
- fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
- fHistDeltaEtaDeltaPhiLeadingJet(0x0),
- fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
- fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
+ fHistEvtClass(0x0),
+ fHistCentrality(0x0),
+ fHistPtJet(new TH1F*[fkNbranches]),
+ fHistEtaPhiJet(new TH2F*[fkNbranches]),
+ fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
+ fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
+ fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
+ fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
+ fHistPtFraction(new TH2F*[fkEvtClasses]),
fHistPtPtExtra(0x0),
fHistPtResponse(new TH2F*[fkEvtClasses]),
fHistPtSmearing(new TH2F*[fkEvtClasses]),
- fHistdR(new TH2F*[fkEvtClasses]),
- fHistArea(new TH2F*[fkEvtClasses])
+ fHistDeltaR(new TH2F*[fkEvtClasses]),
+ fHistArea(new TH2F*[fkEvtClasses]),
+ fHistDeltaArea(new TH2F*[fkEvtClasses])
{
// Constructor
{
// Create histograms
// Called once
+ OpenFile(1);
+ if(!fOutputList) fOutputList = new TList;
+ fOutputList->SetOwner(kTRUE);
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
- fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
+ fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
+
+ fHistEvtClass = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
+ fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
- fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
- Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
- 300, 0., 300.);
- fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
-
- fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
- Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
- 100, -2., 2., 100, 0., 2*TMath::Pi());
+ fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
+ Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
+ 250, 0., 250.);
+ fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
+
+ fHistEtaPhiJet[iBranch] = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
+ Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
+ 48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
+ fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
+ Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
+ 48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
}
- fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
- "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
- 100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());
+
fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
- 300, 0., 300., 300, 0., 300.);
+ 250, 0., 250., 250, 0., 250.);
for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
+ fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
+ "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
+ 101, -1.01, 1.01, 101, -1.01, 1.01);
+ fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
+ "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
+ 101, -1.01, 1.01, 101, -1.01, 1.01);
+ fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
+ "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
+ 100, -2., 2., 100, TMath::Pi(), TMath::Pi());
fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
- 300, 0., 300., 300, 0., 300.);
- fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
- 200, -50., 150., 300, 0., 300.);
- fHistdR[iEvtClass] = new TH2F(Form("hist_dR%i",iEvtClass), "", 200, -50.,150., 200, -1.,1.);
- fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "", 200, -50.,150., 100, 0.,1.);
-
- fHistEtaPhiLeadingJetCut[iEvtClass] = new TH2F(Form("fHistEtaPhiLeadingJetCut%i", iEvtClass),
- Form("#eta - #phi of leading jet from event class %i;#eta;#phi", iEvtClass),
- 100, -2., 2., 100, 0., 2*TMath::Pi());
- fHistDeltaEtaEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaLeadingJet%i", iEvtClass),
- Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#eta", iEvtClass),
- 100, -1., 1., 100, -.5, .5);
- fHistDeltaPtEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaLeadingJet%i", iEvtClass),
- Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#p_{T}", iEvtClass),
- 100, -1., 1., 200, -50., 150);
+ 250, 0., 250., 250, 0., 250.);
+ fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
+ 200, -50., 150., 250, 0., 250.);
+ fHistDeltaR[iEvtClass] = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
+ fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
+ fHistDeltaArea[iEvtClass] = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3);
+
+ fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
+ Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
+ 60, -.6, .6, 100, -.5, .5);
+ fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
+ Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
+ 60, -.6, .6, 200, -50., 150);
+ fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
+ Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
+ 100, 0., 1., 250, 0., 250.);
}
- OpenFile(1);
- fOutputList = new TList;
+
fOutputList->Add(fHistEvtSelection);
+ fOutputList->Add(fHistEvtClass);
+ fOutputList->Add(fHistCentrality);
+
for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
- fOutputList->Add(fHistPtLeadingJet[iBranch]);
- fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
+ fOutputList->Add(fHistPtJet[iBranch]);
+ fOutputList->Add(fHistEtaPhiJet[iBranch]);
+ fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
}
- fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
+
fOutputList->Add(fHistPtPtExtra);
+
for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
+ fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
+ fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
+ fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
+ fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
+ fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
+ fOutputList->Add(fHistPtFraction[iEvtClass]);
+
fOutputList->Add(fHistPtResponse[iEvtClass]);
fOutputList->Add(fHistPtSmearing[iEvtClass]);
- fOutputList->Add(fHistdR[iEvtClass]);
+ fOutputList->Add(fHistDeltaR[iEvtClass]);
fOutputList->Add(fHistArea[iEvtClass]);
- fOutputList->Add(fHistEtaPhiLeadingJetCut[iEvtClass]);
- fOutputList->Add(fHistDeltaEtaEtaLeadingJet[iEvtClass]);
- fOutputList->Add(fHistDeltaPtEtaLeadingJet[iEvtClass]);
+ fOutputList->Add(fHistDeltaArea[iEvtClass]);
}
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
+ TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
+ if (h1){
+ h1->Sumw2();
+ continue;
+ }
+ }
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1, fOutputList);
}
void AliAnalysisTaskJetResponse::UserExec(Option_t *)
return;
}
+ // -- event selection --
fHistEvtSelection->Fill(1); // number of events before event selection
// physics selection
AliInputEventHandler* inputHandler = (AliInputEventHandler*)
((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
- Printf(" Trigger Selection: event REJECTED ... ");
+ if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
fHistEvtSelection->Fill(2);
PostData(1, fOutputList);
return;
if ((nTracksPrim < fMinContribVtx) ||
(primVtx->GetZ() < fVtxZMin) ||
(primVtx->GetZ() > fVtxZMax) ){
- Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
+ if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
fHistEvtSelection->Fill(3);
PostData(1, fOutputList);
return;
// event class selection (from jet helper task)
Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
- Printf("Event class %d", eventClass);
+ if(fDebug) Printf("Event class %d", eventClass);
if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
fHistEvtSelection->Fill(4);
PostData(1, fOutputList);
// centrality selection
AliCentrality *cent = fESD->GetCentrality();
Float_t centValue = cent->GetCentralityPercentile("TRK");
- printf("centrality: %f\n", centValue);
+ if(fDebug) printf("centrality: %f\n", centValue);
if (centValue < fCentMin || centValue > fCentMax){
- fHistEvtSelection->Fill(5);
+ fHistEvtSelection->Fill(4);
PostData(1, fOutputList);
return;
}
fHistEvtSelection->Fill(0); // accepted events
+ fHistEvtClass->Fill(eventClass);
+ fHistCentrality->Fill(centValue);
+ // -- end event selection --
-
+ // fetch jets
TClonesArray *aodJets[2];
aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
- AliAODJet *leadingJet[2] = { 0x0, 0x0 };
for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
fListJets[iJetType]->Clear();
- if (!aodJets[iJetType])
- continue;
+ if (!aodJets[iJetType]) continue;
+
for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
- if (jet)
- fListJets[iJetType]->Add(jet);
+ if (jet) fListJets[iJetType]->Add(jet);
}
fListJets[iJetType]->Sort();
- leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
-
- if (leadingJet[iJetType])
- fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
}
-
- // check if leading jets are close in eta-phi and compare their Pt
- if (leadingJet[0] && leadingJet[1]) {
- Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
- Float_t deltaPhi = TVector2::Phi_mpi_pi(leadingJet[0]->Phi() - leadingJet[1]->Phi());
-
- if (eventClass > -1 && eventClass < fkEvtClasses){
-
- if(leadingJet[1]->Eta()<fJetEtaMax && leadingJet[1]->Eta()>fJetEtaMin){
- fHistEtaPhiLeadingJetCut[eventClass]->Fill(leadingJet[1]->Eta(), leadingJet[1]->Phi());
- }
- }
-
- // leading jets in "jet" acceptance
- if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
- leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
- Printf("Jet not in eta acceptance.");
- }
- else{
- // check association of jets
- fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
-
- if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
- printf("leading jets two far apart\n");
- else {
- Float_t pt0 = leadingJet[0]->Pt();
- Float_t pt1 = leadingJet[1]->Pt();
- Float_t dPt = pt1-pt0;
- Float_t jetarea = leadingJet[1]->EffectiveAreaCharged();
-
- fHistPtLeadingJet[0]->Fill(pt0);
- fHistPtLeadingJet[1]->Fill(pt1);
-
- fHistPtPtExtra->Fill(pt0, pt1);
-
- if (eventClass > -1 && eventClass < fkEvtClasses){
- fHistPtResponse[eventClass]->Fill(pt1, pt0);
- fHistPtSmearing[eventClass]->Fill(dPt, pt0);
-
- Float_t dR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
- fHistdR[eventClass]->Fill(dPt, dR);
- fHistArea[eventClass]->Fill(dPt, jetarea);
-
- fHistDeltaEtaEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), deltaEta);
- fHistDeltaPtEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), dPt);
- }
- }
- }
+
+ // jet matching
+ static TArrayI aMatchIndex(fListJets[0]->GetEntries());
+ static TArrayF aPtFraction(fListJets[0]->GetEntries());
+ if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
+ if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
+
+ // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
+ AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
+ fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
+ aMatchIndex, aPtFraction, fDebug);
+
+ // loop over matched jets
+ Int_t ir = -1; // index of matched reconstruced jet
+ Float_t fraction = 0.;
+ AliAODJet *jet[2] = { 0x0, 0x0 };
+ Float_t jetEta[2] = { 0., 0. };
+ Float_t jetPhi[2] = { 0., 0. };
+ Float_t jetPt[2] = { 0., 0. };
+ Float_t jetArea[2] = { 0., 0. };
+
+ for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
+ ir = aMatchIndex[ig];
+ if(ir<0) continue;
+ fraction = aPtFraction[ig];
+
+ // fetch jets
+ jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
+ jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
+ if(!jet[0] || !jet[1]) continue;
+
+ for(Int_t i=0; i<fkNbranches; ++i){
+ jetEta[i] = jet[i]->Eta();
+ jetPhi[i] = jet[i]->Phi();
+ jetPt[i] = jet[i]->Pt();
+ jetArea[i] = jet[i]->EffectiveAreaCharged();
+ }
+ if(eventClass > -1 && eventClass < fkEvtClasses){
+ fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
+ }
+
+ if(fraction<fJetPtFractionMin) continue;
+
+ // calculate parameters of associated jets
+ Float_t deltaPt = jetPt[1]-jetPt[0];
+ Float_t deltaEta = jetEta[1]-jetEta[0];
+ Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
+ Float_t deltaR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
+ Float_t deltaArea = jetArea[1]-jetArea[0];
+
+ // fill histograms before acceptance cut
+ for(Int_t i=0; i<fkNbranches; ++i){
+ fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
+ }
+ if(eventClass > -1 && eventClass < fkEvtClasses){
+ fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
+ }
+
+ // jet acceptance + minimum pT check
+ if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
+ jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
+
+ if(fDebug){
+ Printf("Jet not in eta acceptance.");
+ Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
+ Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
+ }
+ continue;
+ }
+ if(jetPt[1] < fJetPtMin){
+ if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
+ continue;
+ }
+
+
+
+ // fill histograms
+ for(Int_t i=0; i<fkNbranches; ++i){
+ fHistPtJet[i] -> Fill(jetPt[i]);
+ fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
+ }
+
+ fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
+
+ if(eventClass > -1 && eventClass < fkEvtClasses){
+ fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
+
+ fHistPtResponse[eventClass] -> Fill(jetPt[1], jetPt[0]);
+ fHistPtSmearing[eventClass] -> Fill(deltaPt, jetPt[0]);
+
+ fHistDeltaPtEtaJet[eventClass] -> Fill(jetEta[0], deltaPt);
+ fHistDeltaEtaEtaJet[eventClass] -> Fill(jetEta[0], deltaEta);
+
+ fHistDeltaR[eventClass] -> Fill(deltaPt, deltaR);
+ fHistArea[eventClass] -> Fill(deltaPt, jetArea[1]);
+ fHistDeltaArea[eventClass] -> Fill(deltaPt, deltaArea);
+
+ }
}
PostData(1, fOutputList);
if (!GetOutputData(1))
return;
}
+