#include "AliAnalysisTaskFastEmbedding.h"
#include "AliAnalysisManager.h"
+#include "AliESDtrack.h"
#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
,fAODevent(0)
,fAODtree(0)
,fAODfile(0)
- ,rndm(new TRandom3())
+ ,rndm(0)
,fAODPathArray(0)
,fAODPath("AliAOD.root")
,fTrackBranch("aodExtraTracks")
,fMCparticlesBranch("aodExtraMCparticles")
+ ,fJetBranch("")
,fEntry(0)
- ,fJobId(0)
- ,fNEntriesPerJob(1000)
,fEmbedMode(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)
{
,fAODevent(0)
,fAODtree(0)
,fAODfile(0)
- ,rndm(new TRandom3())
+ ,rndm(0)
,fAODPathArray(0)
,fAODPath("AliAOD.root")
,fTrackBranch("aodExtraTracks")
,fMCparticlesBranch("aodExtraMCparticles")
+ ,fJetBranch("")
,fEntry(0)
- ,fJobId(0)
- ,fNEntriesPerJob(1000)
,fEmbedMode(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());
,fAODPath(copy.fAODPath)
,fTrackBranch(copy.fTrackBranch)
,fMCparticlesBranch(copy.fMCparticlesBranch)
+ ,fJetBranch(copy.fJetBranch)
,fEntry(copy.fEntry)
- ,fJobId(copy.fJobId)
- ,fNEntriesPerJob(copy.fNEntriesPerJob)
,fEmbedMode(copy.fEmbedMode)
,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
}
fAODPath = o.fAODPath;
fTrackBranch = o.fTrackBranch;
fMCparticlesBranch = o.fMCparticlesBranch;
+ fJetBranch = o.fJetBranch;
fEntry = o.fEntry;
- fJobId = o.fJobId;
- fNEntriesPerJob = o.fNEntriesPerJob;
fEmbedMode = o.fEmbedMode;
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()");
-
-
- // embed mode with AOD
- if(fEmbedMode==kAODFull || fEmbedMode==kAODJetTracks || fEmbedMode==kAODJet4Mom){
-
- // open input AOD
- if(fAODPathArray){
- Int_t nFiles = fAODPathArray->GetEntries();
- Int_t n = rndm->Integer(nFiles);
- AliInfo(Form("Select file %d", n));
- TObjString *objStr = (TObjString*) fAODPathArray->At(n);
- if(!objStr){
- AliError("Could not get path of aod file.");
- return;
- }
- fAODPath = objStr->GetString();
- }
-
- TDirectory *owd = gDirectory;
- fAODfile = TFile::Open(fAODPath.Data());
- owd->cd();
- if(!fAODfile){
- AliError("Could not open AOD file.");
- return;
- }
-
- fAODtree = (TTree*)fAODfile->Get("aodTree");
-
- if(!fAODtree){
- AliError("AOD tree not found.");
- return;
- }
- fAODevent = new AliAODEvent();
- fAODevent->ReadFromTree(fAODtree);
- if(!fAODevent){
- AliError("AOD event not found.");
- 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);
}
// Initialization
if(fDebug > 1) Printf("AliAnalysisTaskFastEmbedding::Init()");
-
- // set seed for rndm according to sub-job id *not implemented yet*
}
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;
}
- Int_t maxEntry = fEntry+fNEntriesPerJob-1;
+ // fetch jets
+ TClonesArray *aodJets = 0;
+ 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;
+ }
+
+
Int_t nEvents = fAODtree->GetEntries();
- if(maxEntry>nEvents) maxEntry=nEvents;
Bool_t useEntry = kFALSE;
while(!useEntry){ // protection need, if no event fulfills requierment
- if(fEntry>maxEntry) fEntry=fJobId*fNEntriesPerJob;
+ if(fEntry>=nEvents){
+ fEntry=0;
+ if(!fAODPathArray){
+ AliDebug(AliLog::kDebug, "Last event in AOD reached, start from entry 0 again.");
+ }
+ else {
+ AliDebug(AliLog::kDebug, "Last event in AOD reached, select new AOD file ...");
+
+ 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;
+ }
+ }
+ }
+
fAODtree->GetEvent(fEntry);
// jet pt selection
if(fEvtSelecMode==kEventsJetPt){
- Int_t nJets = fAODevent->GetNJets();
+ Int_t nJets = aodJets->GetEntries();
for(Int_t ij=0; ij<nJets; ++ij){
- AliAODJet *jet = fAODevent->GetJet(ij);
+ AliAODJet *jet = dynamic_cast<AliAODJet*>(aodJets->At(ij));
+ 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];
if(fEmbedMode==kAODJet4Mom){
// loop over jets
- Int_t nJets = fAODevent->GetNJets();
+ Int_t nJets = aodJets->GetEntries();
fh1TrackN->Fill((Float_t)nJets);
for(Int_t ij=0; ij<nJets; ++ij){
- AliAODJet *jet = fAODevent->GetJet(ij);
+ 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);
+ tmpTr->SetFlags(AliESDtrack::kEmbedded);
new((*tracks)[nAODtracks++]) AliAODTrack(*tmpTr);
dummy = (*tracks)[nAODtracks-1];
}
//__________________________________________________________________________
-/* NEEDS TO BE TESTED
Int_t AliAnalysisTaskFastEmbedding::GetJobID()
{
- Int_t id=0;
+ Int_t id=-1;
- const char* env = gSystem->Getenv("ALIENCOUNTER"); // GRID
- if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
+ const char* env = gSystem->Getenv("ALIEN_PROC_ID"); // GRID
+ //if(!env || !strlen(env)) env = gSystem->Getenv("LSB_JOBINDEX"); // GSI
if(env && strlen(env)){
id= atoi(env);
}
else{
AliInfo("Job index not found. Okay if running locally.");
- Int_t nEvents = fAODtree->GetEntries();
- fNEntriesPerJob = nEvents;
- AliInfo(Form("Asuming single job, set entries per job to maximum %d", fNEntriesPerJob));
}
return id;
-}*/
+}
//__________________________________________________________________________
+Int_t AliAnalysisTaskFastEmbedding::SelectAODfile(){
+ Int_t nFiles = fAODPathArray->GetEntries();
+ Int_t n = rndm->Integer(nFiles);
+ AliInfo(Form("Select AOD file %d", n));
+ TObjString *objStr = (TObjString*) fAODPathArray->At(n);
+ if(!objStr){
+ AliError("Could not get path of aod file.");
+ return -1;
+ }
+ fAODPath = objStr->GetString();
+
+ return n;
+}
+
+//__________________________________________________________________________
+
+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->Close();
+ fAODfile = TFile::Open(fAODPath.Data());
+ owd->cd();
+ if(!fAODfile){
+
+ 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");
+
+ if(!fAODtree){
+ AliError("AOD tree not found.");
+ return -1;
+ }
+
+ delete fAODevent;
+ fAODevent = new AliAODEvent();
+ fAODevent->ReadFromTree(fAODtree);
+ return rc; // file position in AOD path array, if array available
+}