#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
+#include "AliAODExtension.h"
#include "AliAODTrack.h"
#include "AliAODJet.h"
#include "AliAODMCParticle.h"
AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
delete fRef;
delete fRandom;
+
+ if(fTCAJetsOut)fTCAJetsOut->Delete();
+ delete fTCAJetsOut;
+ if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
+ delete fTCAJetsOutRan;
+ if(fTCARandomConesOut)fTCARandomConesOut->Delete();
+ delete fTCARandomConesOut;
+ if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
+ delete fTCARandomConesOutRan;
+ if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
+ delete fAODJetBackgroundOut;
}
-AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
+AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster():
+ AliAnalysisTaskSE(),
fAOD(0x0),
fAODExtension(0x0),
fRef(new TRefArray),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
- fUseGlobalSelection(kFALSE),
fUseBackgroundCalc(kFALSE),
+ fEventSelection(kFALSE),
fFilterMask(0),
+ fFilterType(0),
+ fJetTypes(kJet),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
+ fNSkipLeadingCone(0),
fNRandomCones(0),
fAvgTrials(1),
- fExternalWeight(1),
+ fExternalWeight(1),
+ fTrackEtaWindow(0.9),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
- fJetOutputMinPt(1),
+ fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
fJetTriggerPtCut(0),
+ fVtxZCut(8),
+ fVtxR2Cut(1),
fCentCutUp(0),
fCentCutLo(0),
fNonStdBranch(""),
fGhostArea(0.01),
fActiveAreaRepeats(1),
fGhostEtamax(1.5),
+ fTCAJetsOut(0x0),
+ fTCAJetsOutRan(0x0),
+ fTCARandomConesOut(0x0),
+ fTCARandomConesOutRan(0x0),
+ fAODJetBackgroundOut(0x0),
fRandom(0),
fh1Xsec(0x0),
fh1Trials(0x0),
fRef(new TRefArray),
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
- fUseGlobalSelection(kFALSE),
fUseBackgroundCalc(kFALSE),
+ fEventSelection(kFALSE),
fFilterMask(0),
+ fFilterType(0),
+ fJetTypes(kJet),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
+ fNSkipLeadingCone(0),
fNRandomCones(0),
fAvgTrials(1),
fExternalWeight(1),
+ fTrackEtaWindow(0.9),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
- fJetOutputMinPt(1),
+ fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
fJetTriggerPtCut(0),
+ fVtxZCut(8),
+ fVtxR2Cut(1),
fCentCutUp(0),
fCentCutLo(0),
fNonStdBranch(""),
fGhostArea(0.01),
fActiveAreaRepeats(1),
fGhostEtamax(1.5),
+ fTCAJetsOut(0x0),
+ fTCAJetsOutRan(0x0),
+ fTCARandomConesOut(0x0),
+ fTCARandomConesOutRan(0x0),
+ fAODJetBackgroundOut(0x0),
fRandom(0),
fh1Xsec(0x0),
fh1Trials(0x0),
// 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());
-
+
+ if(fJetTypes&kJet){
+ fTCAJetsOut = new TClonesArray("AliAODJet", 0);
+ fTCAJetsOut->SetName(fNonStdBranch.Data());
+ AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
+ }
- TClonesArray *tcaran = new TClonesArray("AliAODJet", 0);
- tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
- AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data());
+ if(fJetTypes&kJetRan){
+ fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
+ fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
+ AddAODBranch("TClonesArray",&fTCAJetsOutRan,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());
+ fAODJetBackgroundOut = new AliAODJetEventBackground();
+ fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
+ AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
}
}
// create the branch for the random cones with the same R
- TString cName = Form("%sRandomCone",fNonStdBranch.Data());
-
+ TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
+
if(fNRandomCones>0){
- if(!AODEvent()->FindListObject(cName.Data())){
- TClonesArray *tcaC = new TClonesArray("AliAODJet", 0);
- tcaC->SetName(cName.Data());
- AddAODBranch("TClonesArray",&tcaC,fNonStdFile.Data());
+ if(fJetTypes&kRC){
+ if(!AODEvent()->FindListObject(cName.Data())){
+ fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
+ fTCARandomConesOut->SetName(cName.Data());
+ AddAODBranch("TClonesArray",&fTCARandomConesOut,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(fJetTypes&kRCRan){
+ cName = Form("%sRandomCone_random",fNonStdBranch.Data());
+ if(!AODEvent()->FindListObject(cName.Data())){
+ fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
+ fTCARandomConesOutRan->SetName(cName.Data());
+ AddAODBranch("TClonesArray",&fTCARandomConesOutRan,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
+ // we identify 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;
- }
- }
+ // case that we have an AOD extension we need can fetch the background maybe from the extended output
+ fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
}
}
if(!fHistList)fHistList = new TList();
fHistList->SetOwner();
+ PostData(1, fHistList); // post data in any case once
Bool_t oldStatus = TH1::AddDirectoryStatus();
TH1::AddDirectory(kFALSE);
//
// Histogram
- const Int_t nBinPt = 200;
+ const Int_t nBinPt = 100;
Double_t binLimitsPt[nBinPt+1];
for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
if(iPt == 0){
binLimitsPt[iPt] = 0.0;
}
else {// 1.0
- binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
+ binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.0;
}
}
}
}
- const int nChMax = 4000;
+ const int nChMax = 5000;
fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
fh1PtJetsLeadingRecInRan = new TH1F("fh1PtJetsLeadingRecInRan","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fh1PtJetConstRec = new TH1F("fh1PtJetsConstRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fh1PtJetConstLeadingRec = new TH1F("fh1PtJetsConstLeadingRec","Rec jets constituents P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
- fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
- fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
- fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
+ fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn",Form("Rec tracks P_T #eta < %1.2f;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
+ fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn",Form("Rec tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),nBinPt,binLimitsPt);
+ fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn",Form("gen tracks P_T #eta < %1.2f ;p_{T} (GeV/c)",fTrackEtaWindow),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);
fHistList->Add(fh1Z);
fHistList->Add(fh1ZSelect);
fHistList->Add(fh1ZPhySel);
- if(fNRandomCones&&fUseBackgroundCalc){
+ if(fNRandomCones>0&&fUseBackgroundCalc){
for(int i = 0;i<3;i++){
fHistList->Add(fh1BiARandomCones[i]);
fHistList->Add(fh1BiARandomConesRan[i]);
fHistList->Add(fh2NConstLeadingPtRan);
fHistList->Add(fh2TracksLeadingJetPhiPtRan);
fHistList->Add(fh2TracksLeadingJetPhiPtWRan);
- }
+ }
// =========== Switch on Sumw2 for all histos ===========
for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
{
- if(fUseGlobalSelection){
- // no selection by the service task, we continue
- if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
- PostData(1, fHistList);
- 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(!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();
- }
- if(fNRandomCones>0){
- 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();
- }
- }
+ if(fTCAJetsOut)fTCAJetsOut->Delete();
+ if(fTCAJetsOutRan)fTCAJetsOutRan->Delete();
+ if(fTCARandomConesOut)fTCARandomConesOut->Delete();
+ if(fTCARandomConesOutRan)fTCARandomConesOutRan->Delete();
+ if(fAODJetBackgroundOut)fAODJetBackgroundOut->Reset();
AliAODJetEventBackground* externalBackground = 0;
if(!externalBackground&&fBackgroundBranch.Length()){
externalBackground = (AliAODJetEventBackground*)(AODEvent()->FindListObject(fBackgroundBranch.Data()));
+ if((!externalBackground)&&fAODExtension)externalBackground = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
}
//
fh1ZPhySel->Fill(zVtx);
}
-
- if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
+ if(fEventSelection){
+ if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
Float_t yvtx = vtxAOD->GetY();
Float_t xvtx = vtxAOD->GetX();
Float_t r2 = yvtx*yvtx+xvtx*xvtx;
- if(TMath::Abs(zVtx)<8.&&r2<1.){ // apply vertex cut later on
+ if(TMath::Abs(zVtx)<fVtxZCut&&r2<fVtxR2Cut){ // apply vertex cut later on
if(physicsSelection){
selectEvent = true;
}
}
- }
- if(fCentCutUp>0){
- if(cent<fCentCutLo||cent>fCentCutUp){
- selectEvent = false;
}
+ if(fCentCutUp>0){
+ if(cent<fCentCutLo||cent>fCentCutUp){
+ selectEvent = false;
+ }
+ }
+ }else{
+ selectEvent = true;
}
-
}
+
+
if(!selectEvent){
PostData(1, fHistList);
return;
// 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 * fRandom->Rndm() - 0.9;
+ Double_t eta = 2.*fTrackEtaWindow * fRandom->Rndm() - fTrackEtaWindow;
Double_t phi = 2.* TMath::Pi() * fRandom->Rndm();
Double_t theta = 2.*TMath::ATan(TMath::Exp(-eta));
}
// fill the tref array, only needed when we write out jets
- if(jarray){
+ if(fTCAJetsOut){
if(i == 0){
fRef->Delete(); // make sure to delete before placement new...
new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
// 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::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::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef);
//range where to compute background
fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax);
+ const vector <fastjet::PseudoJet> &inclusiveJets = clustSeq.inclusive_jets();
+ const vector <fastjet::PseudoJet> &sortedJets = sorted_by_pt(inclusiveJets);
- inclusiveJets = clustSeq.inclusive_jets();
- sortedJets = sorted_by_pt(inclusiveJets);
fh1NJetsRec->Fill(sortedJets.size());
if(externalBackground){
// carefull has to be filled in a task before
// todo, ReArrange to the botom
- pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
+ pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
}
pt = leadingJet.Pt() - pTback;
// correlation of leading jet with tracks
}
-
+ TLorentzVector vecareab;
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();
- if(tmpPt>fJetOutputMinPt&&jarray){// cut on the non-background subtracted...
- aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
+ if(tmpPt>fJetOutputMinPt&&fTCAJetsOut){// cut on the non-background subtracted...
+ aodOutJet = new ((*fTCAJetsOut)[nAodOutJets++]) AliAODJet(tmpRec);
+ aodOutJet->GetRefTracks()->Clear();
Double_t area1 = clustSeq.area(sortedJets[j]);
aodOutJet->SetEffArea(area1,0);
+ fastjet::PseudoJet vecarea=clustSeq.area_4vector(sortedJets[j]);
+ vecareab.SetPxPyPzE(vecarea.px(),vecarea.py(),vecarea.pz(),vecarea.e());
+ aodOutJet->SetVectorAreaCharged(&vecareab);
}
if(tmpPt<0)tmpPt = 0; // avoid negative weights...
fh1PtJetsRecIn->Fill(tmpPt);
- // Fill Spectra with constituents
- vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
+ // Fill Spectra with constituentsemacs
+ const 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(fNRandomCones>0){
- // create a random jet within the acceptance
- Double_t etaMax = 0.8 - 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++){// test for overlap with leading jets
- 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;
- }
- }
- // test for overlap with previous cones to avoid double counting
- for(int iic = 0;iic<ir;iic++){
- AliAODJet *iicone = (AliAODJet*)rConeArray->At(iic);
- if(iicone){
- if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
- 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);
- }
+
+ 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(part->Pt()>fMaxTrackPtInJet)aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
}
- }// if(fUseBackgroundCalc
-
- 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());
- }
-
+ if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
+ }
+
// correlation
Float_t tmpPhi = tmpRec.Phi();
Float_t tmpEta = tmpRec.Eta();
fh2JetsLeadingPhiPtWC[cenClass]->Fill(dPhi,pt,tmpPt);
}
fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
- }
+ }// loop over reconstructed jets
delete recIter;
+
+
+
+ // Add the random cones...
+ if(fNRandomCones>0&&fTCARandomConesOut){
+ // create a random jet within the acceptance
+ Double_t etaMax = fTrackEtaWindow - 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,fNSkipLeadingCone);jj++){// test for overlap with leading jets
+ 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;
+ }
+ }
+ // test for overlap with previous cones to avoid double counting
+ for(int iic = 0;iic<ir;iic++){
+ AliAODJet *iicone = (AliAODJet*)fTCARandomConesOut->At(iic);
+ if(iicone){
+ if(iicone->DeltaR(&tmpRecC)<2.*fRparam){
+ skip = true;
+ break;
+ }
+ }
+ }
+ if(skip)continue;
+ tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
+ if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
+ if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
+ }// loop over random cones creation
+
+
+ // 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(fTCARandomConesOut){
+ for(int ir = 0;ir < fNRandomCones;ir++){
+ AliAODJet *jC = (AliAODJet*)fTCARandomConesOut->At(ir);
+ if(jC&&jC->DeltaR(vp)<fRparam){
+ if(vp->Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+ 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 = 2.*fTrackEtaWindow* fRandom->Rndm() - fTrackEtaWindow;
+ 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(fTCARandomConesOutRan){
+ for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
+ AliAODJet *jC = (AliAODJet*)fTCARandomConesOutRan->At(ir);
+ if(jC&&jC->DeltaR(&vTmpRanR)<fRparam){
+ if(vTmpRanR.Pt()>fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+ jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0);
+ }
+ }
+ }
+ }
+ }// loop over recparticles
+
+ Float_t jetArea = fRparam*fRparam*TMath::Pi();
+ if(fTCARandomConesOut){
+ for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){
+ // rescale the momntum vectors for the random cones
+
+ AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->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.001; // 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(fTCARandomConesOutRan){
+ for(int ir = 0;ir < fTCARandomConesOutRan->GetEntriesFast();ir++){
+ AliAODJet* rC = (AliAODJet*)fTCARandomConesOutRan->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.001;// 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(fNRandomCones
//background estimates:all bckg jets(0) & wo the 2 hardest(1)
- if(evBkg){
+
+
+
+ if(fAODJetBackgroundOut){
vector<fastjet::PseudoJet> jets2=sortedJets;
if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2);
Double_t bkg1=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);
+ clustSeq.get_median_rho_and_sigma(jets2, range, true, bkg1, sigma1, meanarea1, true);
+ fAODJetBackgroundOut->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);
+ fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
// fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
// fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
}
// find the random jets
- vector <fastjet::PseudoJet> inclusiveJetsRan, sortedJetsRan;
+
fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef);
- inclusiveJetsRan = clustSeqRan.inclusive_jets();
- sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
-
// fill the jet information from random track
+ const vector <fastjet::PseudoJet> &inclusiveJetsRan = clustSeqRan.inclusive_jets();
+ const vector <fastjet::PseudoJet> &sortedJetsRan = sorted_by_pt(inclusiveJetsRan);
fh1NJetsRecRan->Fill(sortedJetsRan.size());
Int_t nRecOverRan = inclusiveJetsRan.size();
Int_t nRecRan = inclusiveJetsRan.size();
+
if(inclusiveJetsRan.size()>0){
AliAODJet leadingJet (sortedJetsRan[0].px(), sortedJetsRan[0].py(), sortedJetsRan[0].pz(), sortedJetsRan[0].E());
Float_t pt = leadingJet.Pt();
Int_t iCount = 0;
+ TLorentzVector vecarearanb;
+
for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){
Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i);
while(pt<ptCut&&iCount<nRecRan){
Float_t tmpPt = tmpRec.Pt();
fh1PtJetsRecInRan->Fill(tmpPt);
// Fill Spectra with constituents
- vector<fastjet::PseudoJet> constituents = clustSeqRan.constituents(sortedJetsRan[j]);
+ const vector<fastjet::PseudoJet> &constituents = clustSeqRan.constituents(sortedJetsRan[j]);
fh1NConstRecRan->Fill(constituents.size());
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);
+ if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){
+ aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec);
Double_t arearan=clustSeqRan.area(sortedJetsRan[j]);
-
- aodOutJetRan->SetEffArea(arearan,0); }
-
-
+ aodOutJetRan->GetRefTracks()->Clear();
+ aodOutJetRan->SetEffArea(arearan,0);
+ fastjet::PseudoJet vecarearan=clustSeqRan.area_4vector(sortedJetsRan[j]);
+ vecarearanb.SetPxPyPzE(vecarearan.px(),vecarearan.py(),vecarearan.pz(),vecarearan.e());
+ aodOutJetRan->SetVectorAreaCharged(&vecarearanb);
-
- 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(evBkg){
+ if(fAODJetBackgroundOut){
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);
+ fAODJetBackgroundOut->SetBackground(2,bkg3,sigma3,meanarea3);
// float areaRandomCone = rRandomCone2 *TMath::Pi();
/*
fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone));
*/
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);
+ if(fTCAJetsOut){
+ for(int i = 0;i < fTCAJetsOut->GetEntriesFast();i++){
+ AliAODJet *jet = (AliAODJet*)fTCAJetsOut->At(i);
Float_t ptSub = jet->Pt() - rho *jet->EffectiveAreaCharged();
if(ptSub>=minPt){
select = true;
}
}
}
-
+
if(select){
static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
fh1CentralitySelect->Fill(cent);
aodH->SetFillAOD(kTRUE);
}
}
-
- if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+ if (fDebug > 2){
+ if(fTCAJetsOut)Printf("%s:%d Rec Jets %d",(char*)__FILE__,__LINE__,fTCAJetsOut->GetEntriesFast());
+ if(fTCAJetsOutRan)Printf("%s:%d Rec Jets Ran %d",(char*)__FILE__,__LINE__,fTCAJetsOutRan->GetEntriesFast());
+ if(fTCARandomConesOut)Printf("%s:%d RC %d",(char*)__FILE__,__LINE__,fTCARandomConesOut->GetEntriesFast());
+ if(fTCARandomConesOutRan)Printf("%s:%d RC Ran %d",(char*)__FILE__,__LINE__,fTCARandomConesOutRan->GetEntriesFast());
+ }
PostData(1, fHistList);
}
if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
Int_t iCount = 0;
- if(type==kTrackAOD){
- AliAODEvent *aod = 0;
- if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
- else aod = AODEvent();
- if(!aod){
- return iCount;
+ if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly){
+ if(type!=kTrackAODextraonly) {
+ AliAODEvent *aod = 0;
+ if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ else aod = AODEvent();
+ if(!aod){
+ if(fDebug>2)Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
+ return iCount;
+ }
+ for(int it = 0;it < aod->GetNumberOfTracks();++it){
+ AliAODTrack *tr = aod->GetTrack(it);
+ Bool_t bGood = false;
+ if(fFilterType == 0)bGood = true;
+ else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
+ else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
+ if((fFilterMask>0)&&((!tr->TestFilterBit(fFilterMask)||(!bGood)))){
+ if(fDebug>10)Printf("%s:%d Not matching filter %d/%d %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks(),fFilterMask,tr->GetFilterMap());
+ continue;
+ }
+ if(TMath::Abs(tr->Eta())>fTrackEtaWindow){
+ if(fDebug>10)Printf("%s:%d Not matching eta %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
+ continue;
+ }
+ if(tr->Pt()<fTrackPtCut){
+ if(fDebug>10)Printf("%s:%d Not matching pt %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
+ continue;
+ }
+ if(fDebug>10)Printf("%s:%d MATCHED %d/%d",(char*)__FILE__,__LINE__,it,aod->GetNumberOfTracks());
+ list->Add(tr);
+ iCount++;
+ }
}
- for(int it = 0;it < aod->GetNumberOfTracks();++it){
- AliAODTrack *tr = aod->GetTrack(it);
- if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
- if(TMath::Abs(tr->Eta())>0.9)continue;
- if(tr->Pt()<fTrackPtCut)continue;
- list->Add(tr);
- iCount++;
+ if(type==kTrackAODextra || type==kTrackAODextraonly) {
+ AliAODEvent *aod = 0;
+ if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ else aod = AODEvent();
+
+ if(!aod){
+ return iCount;
+ }
+ TClonesArray *aodExtraTracks = dynamic_cast<TClonesArray*>(aod->FindListObject("aodExtraTracks"));
+ if(!aodExtraTracks)return iCount;
+ for(int it =0; it<aodExtraTracks->GetEntries(); it++) {
+ AliVParticle *track = dynamic_cast<AliVParticle*> ((*aodExtraTracks)[it]);
+ if (!track) continue;
+
+ AliAODTrack *trackAOD = dynamic_cast<AliAODTrack*> (track);
+ if(!trackAOD)continue;
+ Bool_t bGood = false;
+ if(fFilterType == 0)bGood = true;
+ else if(fFilterType == 1)bGood = trackAOD->IsHybridTPCConstrainedGlobal();
+ else if(fFilterType == 2)bGood = trackAOD->IsHybridGlobalConstrainedGlobal();
+ if((fFilterMask>0)&&((!trackAOD->TestFilterBit(fFilterMask)||(!bGood))))continue;
+ if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue;
+ if(trackAOD->Pt()<fTrackPtCut) continue;
+ list->Add(trackAOD);
+ iCount++;
+ }
}
}
else if (type == kTrackKineAll||type == kTrackKineCharged){
list->Add(part);
}
else {
- if(TMath::Abs(part->Eta())>0.9)continue;
+ if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue;
list->Add(part);
}
iCount++;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff