#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
-#include <TRandom.h>
+#include <TRefArray.h>
#include <TFile.h>
#include <TKey.h>
#include <TH1F.h>
#include "AliJetFinder.h"
#include "AliJetHeader.h"
#include "AliJetReader.h"
-#include "AliJetReaderHeader.h"
-#include "AliUA1JetHeaderV1.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliJetKineReaderHeader.h"
#include "AliGenCocktailEventHeader.h"
#include "AliInputEventHandler.h"
-
+#include "AliAODJetEventBackground.h"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
ClassImp(AliAnalysisTaskJetCluster)
+AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
+ delete fRef;
+}
+
AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
fAOD(0x0),
+ fAODExtension(0x0),
+ fRef(new TRefArray),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseGlobalSelection(kFALSE),
+ fUseBackgroundCalc(kFALSE),
fFilterMask(0),
fTrackTypeRec(kTrackUndef),
- fTrackTypeGen(kTrackUndef),
+ fTrackTypeGen(kTrackUndef),
+ fNSkipLeadingRan(0),
fAvgTrials(1),
fExternalWeight(1),
fRecEtaWindow(0.5),
+ fTrackPtCut(0.),
+ fJetOutputMinPt(1),
+ fNonStdBranch(""),
+ fNonStdFile(""),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
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),
fh2TracksLeadingJetPhiPtWRan(0x0),
fHistList(0x0)
{
-
+ for(int i = 0;i<3;i++){
+ fh1BiARandomCones[i] = 0;
+ fh1BiARandomConesRan[i] = 0;
+ }
}
AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
AliAnalysisTaskSE(name),
fAOD(0x0),
+ fAODExtension(0x0),
+ fRef(new TRefArray),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseGlobalSelection(kFALSE),
+ fUseBackgroundCalc(kFALSE),
fFilterMask(0),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
+ fNSkipLeadingRan(0),
fAvgTrials(1),
fExternalWeight(1),
fRecEtaWindow(0.5),
+ fTrackPtCut(0.),
+ fJetOutputMinPt(1),
+ fNonStdBranch(""),
+ fNonStdFile(""),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
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),
fh2TracksLeadingJetPhiPtWRan(0x0),
fHistList(0x0)
{
- DefineOutput(1, TList::Class());
+ for(int i = 0;i<3;i++){
+ fh1BiARandomCones[i] = 0;
+ fh1BiARandomConesRan[i] = 0;
+ }
+ DefineOutput(1, TList::Class());
}
//
+
+
// Connect the AOD
if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
+
+
+ if(fNonStdBranch.Length()!=0)
+ {
+ // only create the output branch if we have a name
+ // Create a new branch for jets...
+ // -> cleared in the UserExec....
+ // here we can also have the case that the brnaches are written to a separate file
+
+ 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(fUseBackgroundCalc){
+ if(!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
+ // we identifay the extension aod event by looking for the branchname
+ AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+ TObjArray* extArray = aodH->GetExtensions();
+ if (extArray) {
+ TIter next(extArray);
+ while ((fAODExtension=(AliAODExtension*)next())){
+ TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
+ if(fDebug>10){
+ Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
+ fAODExtension->GetAOD()->Dump();
+ }
+ if(obj){
+ if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
+ break;
+ }
+ fAODExtension = 0;
+ }
+ }
+ }
+ }
+
+
OpenFile(1);
if(!fHistList)fHistList = new TList();
+ fHistList->SetOwner();
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fh1Nch = new TH1F("fh1Nch","charged multiplicity; N_{ch}",nChMax,-0.5,nChMax-0.5);
+
+ for(int i = 0;i<3;i++){
+ fh1BiARandomCones[i] = new TH1F(Form("fh1BiARandomCones%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
+ fh1BiARandomConesRan[i] = new TH1F(Form("fh1BiARandomConesRan%d",i),";B_{i}^{A} (GeV/c)",200,-100,100);
+ }
+
fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
fh2NRecJetsPtRan = new TH2F("fh2NRecJetsPtRan","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
fHistList->Add(fh1NConstLeadingRecRan);
fHistList->Add(fh1PtJetsRecInRan);
fHistList->Add(fh1Nch);
+ for(int i = 0;i<3;i++){
+ fHistList->Add(fh1BiARandomCones[i]);
+ fHistList->Add(fh1BiARandomConesRan[i]);
+ }
fHistList->Add(fh2NRecJetsPt);
fHistList->Add(fh2NRecTracksPt);
fHistList->Add(fh2NConstPt);
return;
}
+ // 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(fUseBackgroundCalc){
+ 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
//
}
Bool_t selectEvent = false;
- Bool_t physicsSelection = true; // fInputHandler->IsEventSelected();
+ Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
if(fAOD){
const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
- if(vtxAOD->GetNContributors()>0){
+ TString vtxTitle(vtxAOD->GetTitle());
+
+ if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
Float_t zvtx = vtxAOD->GetZ();
Float_t yvtx = vtxAOD->GetY();
Float_t xvtx = vtxAOD->GetX();
vector<fastjet::PseudoJet> inputParticlesRec;
vector<fastjet::PseudoJet> inputParticlesRecRan;
+
+
+ // create a random jet within the acceptance
+ const float rRandomCone2 = 0.4*0.4;
+ float etaRandomCone = gRandom->Rndm()-0.5; // +- 0.5
+ float phiRandomCone = gRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
+ float ptRandomCone = 0;
+ float ptRandomConeRan = 0;
+
for(int i = 0; i < recParticles.GetEntries(); i++){
AliVParticle *vp = (AliVParticle*)recParticles.At(i);
// Carefull energy is not well determined in real data, should not matter for p_T scheme?
jInp.set_user_index(i);
inputParticlesRec.push_back(jInp);
- // the randomized input changes eta and phi, but keeps the p_T
- Double_t pT = vp->Pt();
- Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
- Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
-
-
- Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
- Double_t pZ = pT/TMath::Tan(theta);
-
- Double_t pX = pT * TMath::Cos(phi);
- Double_t pY = pT * TMath::Sin(phi);
- Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
+ if(evBkg){
+ float deta = vp->Eta()-etaRandomCone;
+ float dphi = vp->Phi()-phiRandomCone;
+ if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
+ if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
+ float dr2 = dphi*dphi+deta*deta;
+ if(dr2<rRandomCone2){
+ // within random jet
+ ptRandomCone += vp->Pt();
+ }
+ }
- fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
- jInpRan.set_user_index(i);
- inputParticlesRecRan.push_back(jInpRan);
+ // the randomized input changes eta and phi, but keeps the p_T
+ if(i>=fNSkipLeadingRan){// eventually skip the leading particles
+ Double_t pT = vp->Pt();
+ Double_t eta = 1.8 * gRandom->Rndm() - 0.9;
+ Double_t phi = 2.* TMath::Pi() * gRandom->Rndm();
+
+ Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
+ Double_t pZ = pT/TMath::Tan(theta);
+
+ Double_t pX = pT * TMath::Cos(phi);
+ Double_t pY = pT * TMath::Sin(phi);
+ Double_t p = TMath::Sqrt(pT*pT+pZ*pZ);
+ fastjet::PseudoJet jInpRan(pX,pY,pZ,p);
+ jInpRan.set_user_index(i);
+ inputParticlesRecRan.push_back(jInpRan);
+
+ if(evBkg){
+ float deta = eta-etaRandomCone;
+ float dphi = phi-phiRandomCone;
+ if(dphi>TMath::Pi())dphi -= 2.*TMath::Pi();
+ if(dphi<-TMath::Pi())dphi += 2.*TMath::Pi();
+ float dr2 = dphi*dphi+deta*deta;
+ if(dr2<rRandomCone2){
+ // within random jet
+ ptRandomConeRan += pT;
+ }
+ }
+ }
+ // fill the tref array, only needed when we write out jets
+ if(jarray){
+ if(i == 0){
+ fRef->Delete(); // make sure to delete before placement new...
+ new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
+ }
+ fRef->Add(vp);
+ }
}
if(inputParticlesRec.size()==0){
}
// run fast jet
+ // employ setters for these...
+
+
+ // now create the object that holds info about ghosts
+ if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){
+ // reduce CPU time...
+ fGhostArea = 0.5;
+ fActiveAreaRepeats = 0;
+ }
+ 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);
vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
- fastjet::ClusterSequence clustSeq(inputParticlesRec, jetDef);
+ 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);
// loop over all jets an fill information, first on is the leading jet
- Int_t nRecOver = inclusiveJets.size();
- Int_t nRec = inclusiveJets.size();
- if(inclusiveJets.size()>0){
- AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
+ Int_t nRecOver = inclusiveJets.size();
+ Int_t nRec = inclusiveJets.size();
+ if(inclusiveJets.size()>0){
+ AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
Float_t pt = leadingJet.Pt();
+ Int_t nAodOutJets = 0;
+ Int_t nAodOutTracks = 0;
+ AliAODJet *aodOutJet = 0;
Int_t iCount = 0;
for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
if(phi<0)phi+=TMath::Pi()*2.;
Float_t eta = leadingJet.Eta();
pt = leadingJet.Pt();
-
+
// correlation of leading jet with tracks
TIterator *recIter = recParticles.MakeIterator();
- AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());
-
recIter->Reset();
+ AliVParticle *tmpRecTrack = 0;
while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
Float_t tmpPt = tmpRecTrack->Pt();
// correlation
Float_t dPhi = phi - tmpPhi;
if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
- // Float_t dEta = eta - tmpRecTrack->Eta();
fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
- }
+ }
+
+
+
+ for(int j = 0; j < nRec;j++){
+ AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
+ aodOutJet = 0;
+ nAodOutTracks = 0;
+ Float_t tmpPt = tmpRec.Pt();
+ fh1PtJetsRecIn->Fill(tmpPt);
+ // Fill Spectra with constituents
+ vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
+ fh1NConstRec->Fill(constituents.size());
+ fh2PtNch->Fill(nCh,tmpPt);
+ fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
+ fh2NConstPt->Fill(tmpPt,constituents.size());
+ // loop over constiutents and fill spectrum
+
+ // Add the jet information and the track references to the output AOD
+
+ if(tmpPt>fJetOutputMinPt&&jarray){
+ aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
+ Double_t area=clustSeq.area(sortedJets[j]);
+
+ aodOutJet->SetEffArea(area,0);
+ }
+
+ for(unsigned int ic = 0; ic < constituents.size();ic++){
+ AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
+ fh1PtJetConstRec->Fill(part->Pt());
+ if(aodOutJet){
+ aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
+ }
+ if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
+ }
- for(int j = 0; j < nRec;j++){
- AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
- Float_t tmpPt = tmpRec.Pt();
- fh1PtJetsRecIn->Fill(tmpPt);
- // Fill Spectra with constituents
- vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
- fh1NConstRec->Fill(constituents.size());
- fh2PtNch->Fill(nCh,tmpPt);
- fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
- fh2NConstPt->Fill(tmpPt,constituents.size());
- // loop over constiutents and fill spectrum
- for(unsigned int ic = 0; ic < constituents.size();ic++){
- AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
- fh1PtJetConstRec->Fill(part->Pt());
- if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
- }
- // correlation
- Float_t tmpPhi = tmpRec.Phi();
- Float_t tmpEta = tmpRec.Eta();
- if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
+
+
- if(j==0){
- fh1PtJetsLeadingRecIn->Fill(tmpPt);
- fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
- fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
- fh1NConstLeadingRec->Fill(constituents.size());
- fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
- continue;
- }
- fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
- fh2JetEtaPt->Fill(tmpEta,tmpPt);
- Float_t dPhi = phi - tmpPhi;
- if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
- if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
- Float_t dEta = eta - tmpRec.Eta();
- fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
- fh2JetsLeadingPhiPt->Fill(dPhi,pt);
- fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
- }
- delete recIter;
- }
+
+
+ // correlation
+ Float_t tmpPhi = tmpRec.Phi();
+ Float_t tmpEta = tmpRec.Eta();
+ if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
+
+ if(j==0){
+ fh1PtJetsLeadingRecIn->Fill(tmpPt);
+ fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
+ fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
+ fh1NConstLeadingRec->Fill(constituents.size());
+ fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
+ continue;
+ }
+ fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
+ fh2JetEtaPt->Fill(tmpEta,tmpPt);
+ Float_t dPhi = phi - tmpPhi;
+ if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
+ if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
+ Float_t dEta = eta - tmpRec.Eta();
+ fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
+ fh2JetsLeadingPhiPt->Fill(dPhi,pt);
+ fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
+ }
+ delete recIter;
+
+ //background estimates:all bckg jets(0) & wo the 2 hardest(1)
+ if(evBkg){
+ vector<fastjet::PseudoJet> 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.;
+
+ float areaRandomCone = rRandomCone2 *TMath::Pi();
+ clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
+ evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+ fh1BiARandomCones[0]->Fill(ptRandomCone-(bkg1*areaRandomCone));
+ fh1BiARandomConesRan[0]->Fill(ptRandomConeRan-(bkg1*areaRandomCone));
+
+ clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true);
+ evBkg->SetBackground(1,bkg2,sigma2,meanarea2);
+ fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
+ fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
+ }
+ }
+
- // fill track information
- Int_t nTrackOver = recParticles.GetSize();
+
+
+
+ // fill track information
+ Int_t nTrackOver = recParticles.GetSize();
// do the same for tracks and jets
- if(nTrackOver>0){
+
+ if(nTrackOver>0){
TIterator *recIter = recParticles.MakeIterator();
AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
Float_t pt = tmpRec->Pt();
+
// Printf("Leading track p_t %3.3E",pt);
for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
// find the random jets
vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
- fastjet::ClusterSequence clustSeqRan(inputParticlesRecRan, jetDef);
+ fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
inclusiveJetsRan = clustSeqRan.inclusive_jets();
sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
Float_t dPhi = phi - tmpPhi;
if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
- // Float_t dEta = eta - tmpRecTrack->Eta();
fh2TracksLeadingJetPhiPtRan->Fill(dPhi,pt);
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();
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.;
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, true);
+ evBkg->SetBackground(2,bkg3,sigma3,meanarea3);
+ float areaRandomCone = rRandomCone2 *TMath::Pi();
+ fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
+ fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone));
+ }
+
+
+
}
TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
if(!tca)return iCount;
for(int it = 0;it < tca->GetEntriesFast();++it){
- AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
+ AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
if(!part->IsPhysicalPrimary())continue;
if(type == kTrackAODMCAll){
list->Add(part);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff