#include <TROOT.h>
-#include <TRandom.h>
+#include <TRandom3.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
delete fRef;
+ delete fRandom;
}
AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseGlobalSelection(kFALSE),
+ fUseBackgroundCalc(kFALSE),
fFilterMask(0),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
+ fNRandomCones(10),
fAvgTrials(1),
fExternalWeight(1),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
fJetOutputMinPt(1),
+ fJetTriggerPtCut(0),
+ fCentCutUp(0),
+ fCentCutLo(0),
fNonStdBranch(""),
+ fBackgroundBranch(""),
fNonStdFile(""),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
fGhostArea(0.01),
fActiveAreaRepeats(1),
fGhostEtamax(1.5),
+ fRandom(0),
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
fh1PtJetsRecInRan(0x0),
fh1PtTracksGenIn(0x0),
fh1Nch(0x0),
+ fh1CentralityPhySel(0x0),
+ fh1Centrality(0x0),
+ fh1CentralitySelect(0x0),
+ fh1ZPhySel(0x0),
+ fh1Z(0x0),
+ fh1ZSelect(0x0),
fh2NRecJetsPt(0x0),
fh2NRecTracksPt(0x0),
fh2NConstPt(0x0),
fh2TracksLeadingJetPhiPtWRan(0x0),
fHistList(0x0)
{
-
+ for(int i = 0;i<3;i++){
+ fh1BiARandomCones[i] = 0;
+ fh1BiARandomConesRan[i] = 0;
+ }
+ for(int i = 0;i<kMaxCent;i++){
+ fh2JetsLeadingPhiPtC[i] = 0;
+ fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
+ fh2TracksLeadingJetPhiPtC[i] = 0;
+ fh2TracksLeadingJetPhiPtWC[i] = 0;
+ }
}
AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseGlobalSelection(kFALSE),
+ fUseBackgroundCalc(kFALSE),
fFilterMask(0),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
+ fNRandomCones(5),
fAvgTrials(1),
fExternalWeight(1),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
fJetOutputMinPt(1),
+ fJetTriggerPtCut(0),
+ fCentCutUp(0),
+ fCentCutLo(0),
fNonStdBranch(""),
+ fBackgroundBranch(""),
fNonStdFile(""),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
fGhostArea(0.01),
fActiveAreaRepeats(1),
fGhostEtamax(1.5),
+ fRandom(0),
fh1Xsec(0x0),
fh1Trials(0x0),
fh1PtHard(0x0),
fh1PtJetsRecInRan(0x0),
fh1PtTracksGenIn(0x0),
fh1Nch(0x0),
+ fh1CentralityPhySel(0x0),
+ fh1Centrality(0x0),
+ fh1CentralitySelect(0x0),
+ fh1ZPhySel(0x0),
+ fh1Z(0x0),
+ fh1ZSelect(0x0),
fh2NRecJetsPt(0x0),
fh2NRecTracksPt(0x0),
fh2NConstPt(0x0),
fh2TracksLeadingJetPhiPtWRan(0x0),
fHistList(0x0)
{
- DefineOutput(1, TList::Class());
+ for(int i = 0;i<3;i++){
+ fh1BiARandomCones[i] = 0;
+ fh1BiARandomConesRan[i] = 0;
+ }
+ for(int i = 0;i<kMaxCent;i++){
+ fh2JetsLeadingPhiPtC[i] = 0;
+ fh2JetsLeadingPhiPtWC[i] = 0; //! jet correlation with leading jet
+ fh2TracksLeadingJetPhiPtC[i] = 0;
+ fh2TracksLeadingJetPhiPtWC[i] = 0;
+ }
+ DefineOutput(1, TList::Class());
}
// Create the output container
//
-
+ fRandom = new TRandom3(0);
// Connect the AOD
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(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());
+ }
+ // create the branch for the random cones with the same R
+ TString cName = Form("%sRandomCone",fNonStdBranch.Data());
+ if(!AODEvent()->FindListObject(cName.Data())){
+ TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
+ tcaC->SetName(cName.Data());
+ AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
+ }
+
+ // create the branch with the random for the random cones on the random event
+ cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+ if(!AODEvent()->FindListObject(cName.Data())){
+ TClonesArray *tcaCran = new TClonesArray("AliAODJet", 0);
+ tcaCran->SetName(cName.Data());
+ AddAODBranch("TClonesArray",&tcaCran,fNonStdFile.Data());
+ }
+ }
if(fNonStdFile.Length()!=0){
//
}
- OpenFile(1);
if(!fHistList)fHistList = new TList();
+ fHistList->SetOwner();
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
binLimitsPt[iPt] = 0.0;
}
else {// 1.0
- binLimitsPt[iPt] = binLimitsPt[iPt-1] + 0.5;
+ binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
}
}
}
}
- const int nChMax = 100;
+ const int nChMax = 4000;
fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
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);
+ fh1Centrality = new TH1F("fh1Centrality",";cent (%)",111,-0.5,110.5);
+ fh1CentralitySelect = new TH1F("fh1CentralitySelect",";cent (%)",111,-0.5,110.5);
+ fh1CentralityPhySel = new TH1F("fh1CentralityPhySel",";cent (%)",111,-0.5,110.5);
+
+ fh1Z = new TH1F("fh1Z",";zvtx",100,-25,25);
+ fh1ZSelect = new TH1F("fh1ZSelect",";zvtx",100,-25,25);
+ fh1ZPhySel = new TH1F("fh1ZPhySel",";zvtx",100,-25,25);
+
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);
nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+ if(fUseBackgroundCalc){
+ 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);
+ }
+ }
+ for(int i = 0;i < kMaxCent;i++){
+ fh2JetsLeadingPhiPtC[i] = (TH2F*)fh2JetsLeadingPhiPt->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPt->GetName(),i+1));
+ fh2JetsLeadingPhiPtWC[i]= (TH2F*)fh2JetsLeadingPhiPtW->Clone(Form("%s_C%02d",fh2JetsLeadingPhiPtW->GetName(),i+1));
+ fh2TracksLeadingJetPhiPtC[i] = (TH2F*)fh2TracksLeadingJetPhiPt->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPt->GetName(),i+1));
+ fh2TracksLeadingJetPhiPtWC[i] = (TH2F*)fh2TracksLeadingJetPhiPtW->Clone(Form("%s_C%02d",fh2TracksLeadingJetPhiPtW->GetName(),i+1));
+ }
const Int_t saveLevel = 3; // large save level more histos
if(saveLevel>0){
fHistList->Add(fh1NConstLeadingRecRan);
fHistList->Add(fh1PtJetsRecInRan);
fHistList->Add(fh1Nch);
+ fHistList->Add(fh1Centrality);
+ fHistList->Add(fh1CentralitySelect);
+ fHistList->Add(fh1CentralityPhySel);
+ fHistList->Add(fh1Z);
+ fHistList->Add(fh1ZSelect);
+ fHistList->Add(fh1ZPhySel);
+ if(fUseBackgroundCalc){
+ for(int i = 0;i<3;i++){
+ fHistList->Add(fh1BiARandomCones[i]);
+ fHistList->Add(fh1BiARandomConesRan[i]);
+ }
+ }
+ for(int i = 0;i < kMaxCent;i++){
+ fHistList->Add(fh2JetsLeadingPhiPtC[i]);
+ fHistList->Add(fh2JetsLeadingPhiPtWC[i]);
+ fHistList->Add(fh2TracksLeadingJetPhiPtC[i]);
+ fHistList->Add(fh2TracksLeadingJetPhiPtWC[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;
+ TClonesArray* rConeArray = 0;
+ TClonesArray* rConeArrayRan = 0;
AliAODJetEventBackground* evBkg = 0;
if(fNonStdBranch.Length()!=0) {
if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
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();
-
+
+ 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();
+
+ TString cName = Form("%sRandomCone",fNonStdBranch.Data());
+ if(AODEvent())rConeArray = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
+ if(!rConeArray)rConeArray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
+ if(rConeArray)rConeArray->Delete();
+
+ cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+ if(AODEvent())rConeArrayRan = (TClonesArray*)(AODEvent()->FindListObject(cName.Data()));
+ if(!rConeArrayRan)rConeArrayRan = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(cName.Data()));
+ if(rConeArrayRan)rConeArrayRan->Delete();
+ }
}
-
-
+ static AliAODJetEventBackground* externalBackground = 0;
+ if(!externalBackground&&fBackgroundBranch.Length()){
+ externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
+ Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
+ }
//
// Execute analysis for current event
//
Bool_t selectEvent = false;
Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
+ Float_t cent = 0;
+ Float_t zVtx = 0;
+ Int_t cenClass = -1;
if(fAOD){
const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
TString vtxTitle(vtxAOD->GetTitle());
+ zVtx = vtxAOD->GetZ();
+
+ cent = fAOD->GetHeader()->GetCentrality();
+ if(cent<10)cenClass = 0;
+ else if(cent<30)cenClass = 1;
+ else if(cent<50)cenClass = 2;
+ else if(cent<80)cenClass = 3;
+ if(physicsSelection){
+ fh1CentralityPhySel->Fill(cent);
+ fh1ZPhySel->Fill(zVtx);
+ }
+
if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
- Float_t zvtx = vtxAOD->GetZ();
Float_t yvtx = vtxAOD->GetY();
Float_t xvtx = vtxAOD->GetX();
Float_t r2 = yvtx*yvtx+xvtx*xvtx;
- if(TMath::Abs(zvtx)<8.&&r2<1.){
- if(physicsSelection)selectEvent = true;
+ if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
+ if(physicsSelection){
+ selectEvent = true;
+ }
}
}
+ if(fCentCutUp>0){
+ if(cent<fCentCutLo||cent>fCentCutUp){
+ selectEvent = false;
+ }
+ }
+
}
if(!selectEvent){
PostData(1, fHistList);
return;
}
-
+ fh1Centrality->Fill(cent);
+ fh1Z->Fill(zVtx);
fh1Trials->Fill("#sum{ntrials}",1);
vector<fastjet::PseudoJet> inputParticlesRec;
vector<fastjet::PseudoJet> inputParticlesRecRan;
+
+ // Generate the random cones
+
+ AliAODJet vTmpRan(1,0,0,1);
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?
- // we talk total momentum here
+ // we take total momentum here
fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
jInp.set_user_index(i);
inputParticlesRec.push_back(jInp);
// 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 eta = 1.8 * fRandom->Rndm() - 0.9;
+ Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
- Double_t theta = 2.*TMath::ATan(TMath::Exp(-2.*eta));
+ Double_t theta = 2.*TMath::ATan(TMath::Exp(-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);
+ vTmpRan.SetPxPyPzE(pX,pY,pZ,p);
}
// fill the tref array, only needed when we write out jets
}
fRef->Add(vp);
}
- }
+ }// recparticles
if(inputParticlesRec.size()==0){
if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__);
// employ setters for these...
- // now create the object that holds info about ghosts
- fastjet::GhostedAreaSpec ghostSpec(fGhostEtamax, fActiveAreaRepeats, fGhostArea);
+ // 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);
+ fastjet::JetDefinition jetDef(fAlgorithm, fRparam, fRecombScheme, fStrategy);
vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
Int_t nRec = inclusiveJets.size();
if(inclusiveJets.size()>0){
AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
+ Double_t area = clustSeq.area(sortedJets[0]);
+ leadingJet.SetEffArea(area,0);
Float_t pt = leadingJet.Pt();
Int_t nAodOutJets = 0;
Int_t nAodOutTracks = 0;
Float_t phi = leadingJet.Phi();
if(phi<0)phi+=TMath::Pi()*2.;
Float_t eta = leadingJet.Eta();
- pt = leadingJet.Pt();
-
+ Float_t pTback = 0;
+ if(externalBackground){
+ // carefull has to be filled in a task before
+ // todo, ReArrange to the botom
+ pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
+ }
+ pt = leadingJet.Pt() - pTback;
// correlation of leading jet with tracks
TIterator *recIter = recParticles.MakeIterator();
recIter->Reset();
if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
+ if(cenClass>=0){
+ fh2TracksLeadingJetPhiPtC[cenClass]->Fill(dPhi,pt);
+ fh2TracksLeadingJetPhiPtWC[cenClass]->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(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();
+ Float_t tmpPtBack = 0;
+ if(externalBackground){
+ // carefull has to be filled in a task before
+ // todo, ReArrange to the botom
+ tmpPtBack = externalBackground->GetBackground(2)*tmpRec.EffectiveAreaCharged();
+ }
+ tmpPt = tmpPt - tmpPtBack;
+ if(tmpPt<0)tmpPt = 0; // avoid negative weights...
+
+ 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 area1 = clustSeq.area(sortedJets[j]);
+ aodOutJet->SetEffArea(area1,0);
+ }
+
+
+ if(fUseBackgroundCalc){
+
+ // create a random jet within the acceptance
+ Double_t etaMax = 0.9 - fRparam;
+ Int_t nCone = 0;
+ Int_t nConeRan = 0;
+ Double_t pTC = 1; // small number
+ for(int ir = 0;ir < fNRandomCones;ir++){
+ Double_t etaC = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax
+ Double_t phiC = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi
+ // massless jet
+ Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
+ Double_t pZC = pTC/TMath::Tan(thetaC);
+ Double_t pXC = pTC * TMath::Cos(phiC);
+ Double_t pYC = pTC * TMath::Sin(phiC);
+ Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
+ AliAODJet tmpRecC (pXC,pYC,pZC, pC);
+ bool skip = false;
+ for(int jj = 0; jj < TMath::Min(nRec,2);jj++){
+ AliAODJet jet (sortedJets[jj].px(), sortedJets[jj].py(), sortedJets[jj].pz(), sortedJets[jj].E());
+ if(jet.DeltaR(& tmpRecC)<2.*fRparam+0.2){
+ skip = true;
+ break;
+ }
+ }
+ if(skip)continue;
+ tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
+ if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRecC);
+ if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRecC);
+ }// random cones
+
+
+ // loop over the reconstructed particles and add up the pT in the random cones
+ // maybe better to loop over randomized particles not in the real jets...
+ // but this by definition brings dow average energy in the whole event
+ AliAODJet vTmpRanR(1,0,0,1);
+ for(int i = 0; i < recParticles.GetEntries(); i++){
+ AliVParticle *vp = (AliVParticle*)recParticles.At(i);
+ if(rConeArray){
+ for(int ir = 0;ir < fNRandomCones;ir++){
+ AliAODJet *jC = (AliAODJet*)rConeArray->At(ir);
+ if(jC&&jC->DeltaR(vp)<fRparam){
+ jC->SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0);
+ }
+ }
+ }// add up energy in cone
+
+ // the randomized input changes eta and phi, but keeps the p_T
+ if(i>=fNSkipLeadingRan){// eventually skip the leading particles
+ Double_t pTR = vp->Pt();
+ Double_t etaR = 1.8 * fRandom->Rndm() - 0.9;
+ Double_t phiR = 2.* TMath::Pi() * fRandom->Rndm();
+
+ Double_t thetaR = 2.*TMath::ATan(TMath::Exp(-etaR));
+ Double_t pZR = pTR/TMath::Tan(thetaR);
+
+ Double_t pXR = pTR * TMath::Cos(phiR);
+ Double_t pYR = pTR * TMath::Sin(phiR);
+ Double_t pR = TMath::Sqrt(pTR*pTR+pZR*pZR);
+ vTmpRanR.SetPxPyPzE(pXR,pYR,pZR,pR);
+ if(rConeArrayRan){
+ for(int ir = 0;ir < fNRandomCones;ir++){
+ AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir);
+ if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
+ jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
+ }
+ }
+ }
+ }
+ }// loop over recparticles
+
+ Float_t jetArea = fRparam*fRparam*TMath::Pi();
+ for(int ir = 0;ir < fNRandomCones;ir++){
+ // rescale the momntum vectors for the random cones
+ if(!rConeArray)continue;
+ AliAODJet *rC = (AliAODJet*)rConeArray->At(ir);
+ if(rC){
+ Double_t etaC = rC->Eta();
+ Double_t phiC = rC->Phi();
+ // massless jet, unit vector
+ pTC = rC->ChargedBgEnergy();
+ if(pTC<=0)pTC = 0.1; // for almost empty events
+ Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
+ Double_t pZC = pTC/TMath::Tan(thetaC);
+ Double_t pXC = pTC * TMath::Cos(phiC);
+ Double_t pYC = pTC * TMath::Sin(phiC);
+ Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
+ rC->SetPxPyPzE(pXC,pYC,pZC, pC);
+ rC->SetBgEnergy(0,0);
+ rC->SetEffArea(jetArea,0);
+ }
+ rC = (AliAODJet*)rConeArrayRan->At(ir);
+ // same wit random
+ if(rC){
+ Double_t etaC = rC->Eta();
+ Double_t phiC = rC->Phi();
+ // massless jet, unit vector
+ pTC = rC->ChargedBgEnergy();
+ if(pTC<=0)pTC = 0.1;// for almost empty events
+ Double_t thetaC = 2.*TMath::ATan(TMath::Exp(-etaC));
+ Double_t pZC = pTC/TMath::Tan(thetaC);
+ Double_t pXC = pTC * TMath::Cos(phiC);
+ Double_t pYC = pTC * TMath::Sin(phiC);
+ Double_t pC = TMath::Sqrt(pTC*pTC+pZC*pZC);
+ rC->SetPxPyPzE(pXC,pYC,pZC, pC);
+ rC->SetBgEnergy(0,0);
+ rC->SetEffArea(jetArea,0);
+ }
+ }
+ }// if(fUseBackgroundCalc
for(unsigned int ic = 0; ic < constituents.size();ic++){
AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
}
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);
Float_t dEta = eta - tmpRec.Eta();
fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
fh2JetsLeadingPhiPt->Fill(dPhi,pt);
+ if(cenClass>=0){
+ fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt);
+ fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
+ }
fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
- }
+ }
delete recIter;
-
- //background estimates:all bckg jets(0) & wo the 2 hardest(1)
+
+ //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 bkg2=0;
Double_t sigma2=0.;
Double_t meanarea2=0.;
-
+
clustSeq.get_median_rho_and_sigma(sortedJets, range, false, bkg1, sigma1, meanarea1, true);
evBkg->SetBackground(0,bkg1,sigma1,meanarea1);
+
+ // fh1BiARandomCones[0]->Fill(omCone-(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);
// fh1PtJetConstRec->Fill(part->Pt());
if(aodOutJetRan){
aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index()));
- }}
+ }
+ }
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));
+ */
}
}
-
+
+
+ // do the event selection if activated
+ if(fJetTriggerPtCut>0){
+ bool select = false;
+ Float_t minPt = fJetTriggerPtCut;
+ /*
+ // hard coded for now ...
+ // 54.50 44.5 29.5 18.5 for anti-kt rejection 1E-3
+ if(cent<10)minPt = 50;
+ else if(cent<30)minPt = 42;
+ else if(cent<50)minPt = 28;
+ else if(cent<80)minPt = 18;
+ */
+ float rho = 0;
+ if(externalBackground)rho = externalBackground->GetBackground(2);
+ if(jarray){
+ for(int i = 0;i < jarray->GetEntriesFast();i++){
+ AliAODJet *jet = (AliAODJet*)jarray->At(i);
+ Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
+ if(ptSub>=minPt){
+ select = true;
+ break;
+ }
+ }
+ }
+
+ if(select){
+ static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+ fh1CentralitySelect->Fill(cent);
+ fh1ZSelect->Fill(zVtx);
+ aodH->SetFillAOD(kTRUE);
+ }
+ }
if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
PostData(1, fHistList);
AliAODTrack *tr = aod->GetTrack(it);
if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
if(TMath::Abs(tr->Eta())>0.9)continue;
- // if(tr->Pt()<0.3)continue;
+ if(tr->Pt()<fTrackPtCut)continue;
list->Add(tr);
iCount++;
}
if(!mcEvent->IsPhysicalPrimary(it))continue;
AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
if(type == kTrackKineAll){
+ if(part->Pt()<fTrackPtCut)continue;
list->Add(part);
iCount++;
}
else if(type == kTrackKineCharged){
if(part->Particle()->GetPDG()->Charge()==0)continue;
+ if(part->Pt()<fTrackPtCut)continue;
list->Add(part);
iCount++;
}
AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it));
if(!part->IsPhysicalPrimary())continue;
if(type == kTrackAODMCAll){
+ if(part->Pt()<fTrackPtCut)continue;
list->Add(part);
iCount++;
}
else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
if(part->Charge()==0)continue;
+ if(part->Pt()<fTrackPtCut)continue;
if(kTrackAODMCCharged){
list->Add(part);
}
}// AODMCparticle
list->Sort();
return iCount;
-
}
/*
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff