#include <TH2F.h>
#include <TH3F.h>
#include <TProfile.h>
+#include <TF1.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
// get info on how fastjet was configured
#include "fastjet/config.h"
+using std::vector;
ClassImp(AliAnalysisTaskJetCluster)
AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
+ //
+ // Destructor
+ //
+
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;
}
fUseBackgroundCalc(kFALSE),
fEventSelection(kFALSE),
fFilterMask(0),
+ fFilterMaskBestPt(0),
fFilterType(0),
+ fJetTypes(kJet),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
fJetTriggerPtCut(0),
fVtxZCut(8),
fVtxR2Cut(1),
fCentCutUp(0),
fCentCutLo(0),
+ fStoreRhoLeadingTrackCorr(kFALSE),
fNonStdBranch(""),
fBackgroundBranch(""),
fNonStdFile(""),
+ fMomResH1(0x0),
+ fMomResH2(0x0),
+ fMomResH3(0x0),
+ fMomResH1Fit(0x0),
+ fMomResH2Fit(0x0),
+ fMomResH3Fit(0x0),
+ fhEffH1(0x0),
+ fhEffH2(0x0),
+ fhEffH3(0x0),
+ fUseTrPtResolutionSmearing(kFALSE),
+ fUseDiceEfficiency(kFALSE),
+ fDiceEfficiencyMinPt(0.),
+ fUseTrPtResolutionFromOADB(kFALSE),
+ fUseTrEfficiencyFromOADB(kFALSE),
+ fPathTrPtResolution(""),
+ fPathTrEfficiency(""),
+ fChangeEfficiencyFraction(0.),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
fStrategy(fastjet::Best),
fh2PtNchNRan(0x0),
fh2TracksLeadingJetPhiPtRan(0x0),
fh2TracksLeadingJetPhiPtWRan(0x0),
+ fh2CentvsRho(0x0),
+ fh2CentvsSigma(0x0),
+ fh2MultvsRho(0x0),
+ fh2MultvsSigma(0x0),
+ fh3CentvsRhoLeadingTrackPtQ1(0x0),
+ fh3CentvsRhoLeadingTrackPtQ2(0x0),
+ fh3CentvsRhoLeadingTrackPtQ3(0x0),
+ fh3CentvsRhoLeadingTrackPtQ4(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ1(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ2(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ3(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ4(0x0),
+ fh3MultvsRhoLeadingTrackPtQ1(0x0),
+ fh3MultvsRhoLeadingTrackPtQ2(0x0),
+ fh3MultvsRhoLeadingTrackPtQ3(0x0),
+ fh3MultvsRhoLeadingTrackPtQ4(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ1(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ2(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ3(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ4(0x0),
+ fh2PtGenPtSmeared(0x0),
+ fp1Efficiency(0x0),
+ fp1PtResolution(0x0),
fHistList(0x0)
{
+ //
+ // Constructor
+ //
+
for(int i = 0;i<3;i++){
fh1BiARandomCones[i] = 0;
fh1BiARandomConesRan[i] = 0;
fUseAODTrackInput(kFALSE),
fUseAODMCInput(kFALSE),
fUseBackgroundCalc(kFALSE),
- fEventSelection(kFALSE),
- fFilterMask(0),
+ fEventSelection(kFALSE), fFilterMask(0),
+ fFilterMaskBestPt(0),
fFilterType(0),
+ fJetTypes(kJet),
fTrackTypeRec(kTrackUndef),
fTrackTypeGen(kTrackUndef),
fNSkipLeadingRan(0),
fRecEtaWindow(0.5),
fTrackPtCut(0.),
fJetOutputMinPt(0.150),
+ fMaxTrackPtInJet(100.),
fJetTriggerPtCut(0),
fVtxZCut(8),
fVtxR2Cut(1),
fCentCutUp(0),
fCentCutLo(0),
+ fStoreRhoLeadingTrackCorr(kFALSE),
fNonStdBranch(""),
fBackgroundBranch(""),
fNonStdFile(""),
+ fMomResH1(0x0),
+ fMomResH2(0x0),
+ fMomResH3(0x0),
+ fMomResH1Fit(0x0),
+ fMomResH2Fit(0x0),
+ fMomResH3Fit(0x0),
+ fhEffH1(0x0),
+ fhEffH2(0x0),
+ fhEffH3(0x0),
+ fUseTrPtResolutionSmearing(kFALSE),
+ fUseDiceEfficiency(kFALSE),
+ fDiceEfficiencyMinPt(0.),
+ fUseTrPtResolutionFromOADB(kFALSE),
+ fUseTrEfficiencyFromOADB(kFALSE),
+ fPathTrPtResolution(""),
+ fPathTrEfficiency(""),
+ fChangeEfficiencyFraction(0.),
fRparam(1.0),
fAlgorithm(fastjet::kt_algorithm),
fStrategy(fastjet::Best),
fh2PtNchNRan(0x0),
fh2TracksLeadingJetPhiPtRan(0x0),
fh2TracksLeadingJetPhiPtWRan(0x0),
+ fh2CentvsRho(0x0),
+ fh2CentvsSigma(0x0),
+ fh2MultvsRho(0x0),
+ fh2MultvsSigma(0x0),
+ fh3CentvsRhoLeadingTrackPtQ1(0x0),
+ fh3CentvsRhoLeadingTrackPtQ2(0x0),
+ fh3CentvsRhoLeadingTrackPtQ3(0x0),
+ fh3CentvsRhoLeadingTrackPtQ4(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ1(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ2(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ3(0x0),
+ fh3CentvsSigmaLeadingTrackPtQ4(0x0),
+ fh3MultvsRhoLeadingTrackPtQ1(0x0),
+ fh3MultvsRhoLeadingTrackPtQ2(0x0),
+ fh3MultvsRhoLeadingTrackPtQ3(0x0),
+ fh3MultvsRhoLeadingTrackPtQ4(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ1(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ2(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ3(0x0),
+ fh3MultvsSigmaLeadingTrackPtQ4(0x0),
+ fh2PtGenPtSmeared(0x0),
+ fp1Efficiency(0x0),
+ fp1PtResolution(0x0),
fHistList(0x0)
{
+ //
+ // named ctor
+ //
+
for(int i = 0;i<3;i++){
fh1BiARandomCones[i] = 0;
fh1BiARandomConesRan[i] = 0;
// -> cleared in the UserExec....
// here we can also have the case that the brnaches are written to a separate file
- fTCAJetsOut = new TClonesArray("AliAODJet", 0);
- fTCAJetsOut->SetName(fNonStdBranch.Data());
- AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
-
-
-
- fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
- fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
- AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
-
+ if(fJetTypes&kJet){
+ fTCAJetsOut = new TClonesArray("AliAODJet", 0);
+ fTCAJetsOut->SetName(fNonStdBranch.Data());
+ AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data());
+ }
+
+ if(fJetTypes&kJetRan){
+ fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
+ fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
+ if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
+ fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dFr%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
+ AddAODBranch("TClonesArray",&fTCAJetsOutRan,fNonStdFile.Data());
+ }
+
if(fUseBackgroundCalc){
if(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
fAODJetBackgroundOut = new AliAODJetEventBackground();
fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
+ if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
+ fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dFr%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.)));
+
AddAODBranch("AliAODJetEventBackground",&fAODJetBackgroundOut,fNonStdFile.Data());
}
}
// create the branch for the random cones with the same R
TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone);
-
+ if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
+ cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
+
if(fNRandomCones>0){
- if(!AODEvent()->FindListObject(cName.Data())){
- fTCARandomConesOut = new TClonesArray("AliAODJet", 0);
- fTCARandomConesOut->SetName(cName.Data());
- AddAODBranch("TClonesArray",&fTCARandomConesOut,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())){
- fTCARandomConesOutRan = new TClonesArray("AliAODJet", 0);
- fTCARandomConesOutRan->SetName(cName.Data());
- AddAODBranch("TClonesArray",&fTCARandomConesOutRan,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());
+ }
}
}
fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+ if(fStoreRhoLeadingTrackCorr) {
+ fh2CentvsRho = new TH2F("fh2CentvsRho","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.);
+ fh2CentvsSigma = new TH2F("fh2CentvsSigma","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.);
+ fh2MultvsRho = new TH2F("fh2MultvsRho","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.);
+ fh2MultvsSigma = new TH2F("fh2MultvsSigma","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.);
+
+
+ fh3CentvsRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
+ fh3CentvsRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
+ fh3CentvsRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
+ fh3CentvsRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho; leading p_{t}^{track}", 50,0.,100.,500,0.,250.,100,0.,100.);
+
+ fh3CentvsSigmaLeadingTrackPtQ1 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ1","centrality vs background density Q1; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
+ fh3CentvsSigmaLeadingTrackPtQ2 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ2","centrality vs background density Q2; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
+ fh3CentvsSigmaLeadingTrackPtQ3 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ3","centrality vs background density Q3; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
+ fh3CentvsSigmaLeadingTrackPtQ4 = new TH3F("fh3CentvsSigmaLeadingTrackPtQ4","centrality vs background density Q4; centrality; #sigma(#rho); leading p_{t}^{track}", 50,0.,100.,50,0.,50.,100,0.,100.);
+
+ fh3MultvsRhoLeadingTrackPtQ1 = new TH3F("fh3MultvsRhoLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
+ fh3MultvsRhoLeadingTrackPtQ2 = new TH3F("fh3MultvsRhoLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
+ fh3MultvsRhoLeadingTrackPtQ3 = new TH3F("fh3MultvsRhoLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
+ fh3MultvsRhoLeadingTrackPtQ4 = new TH3F("fh3MultvsRhoLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #rho; leading p_{t}^{track}", 100,0.,5000.,500,0.,250.,100,0.,100.);
+
+ fh3MultvsSigmaLeadingTrackPtQ1 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ1","multiplicity vs background density Q1; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
+ fh3MultvsSigmaLeadingTrackPtQ2 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ2","multiplicity vs background density Q2; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
+ fh3MultvsSigmaLeadingTrackPtQ3 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ3","multiplicity vs background density Q3; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
+ fh3MultvsSigmaLeadingTrackPtQ4 = new TH3F("fh3MultvsSigmaLeadingTrackPtQ4","multiplicity vs background density Q4; multiplicity; #sigma(#rho); leading p_{t}^{track}", 100,0.,5000.,50,0.,50.,100,0.,100.);
+
+ fHistList->Add(fh2CentvsRho);
+ fHistList->Add(fh2CentvsSigma);
+ fHistList->Add(fh2MultvsRho);
+ fHistList->Add(fh2MultvsSigma);
+
+ fHistList->Add(fh3CentvsRhoLeadingTrackPtQ1);
+ fHistList->Add(fh3CentvsRhoLeadingTrackPtQ2);
+ fHistList->Add(fh3CentvsRhoLeadingTrackPtQ3);
+ fHistList->Add(fh3CentvsRhoLeadingTrackPtQ4);
+
+ fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ1);
+ fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ2);
+ fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ3);
+ fHistList->Add(fh3CentvsSigmaLeadingTrackPtQ4);
+
+ fHistList->Add(fh3MultvsRhoLeadingTrackPtQ1);
+ fHistList->Add(fh3MultvsRhoLeadingTrackPtQ2);
+ fHistList->Add(fh3MultvsRhoLeadingTrackPtQ3);
+ fHistList->Add(fh3MultvsRhoLeadingTrackPtQ4);
+
+ fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ1);
+ fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ2);
+ fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ3);
+ fHistList->Add(fh3MultvsSigmaLeadingTrackPtQ4);
+
+ }
+
+ //Detector level effects histos
+ fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+
+ fp1Efficiency = new TProfile("fp1Efficiency","fp1Efficiency",nBinPt,binLimitsPt);
+ fp1PtResolution = new TProfile("fp1PtResolution","fp1PtResolution",nBinPt,binLimitsPt);
+
+ fHistList->Add(fh2PtGenPtSmeared);
+ fHistList->Add(fp1Efficiency);
+ fHistList->Add(fp1PtResolution);
if(fNRandomCones>0&&fUseBackgroundCalc){
for(int i = 0;i<3;i++){
TH1::AddDirectory(oldStatus);
}
-void AliAnalysisTaskJetCluster::Init()
+void AliAnalysisTaskJetCluster::LocalInit()
{
//
// Initialization
if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
+ if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
+ if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
+
+
+ FitMomentumResolution();
+
}
void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODTrackInput);
return;
}
- // fethc the header
+ // fetch the header
}
else{
// assume that the AOD is in the general output...
fESD = dynamic_cast<AliESDEvent*> (InputEvent());
}
}
+
+ //Check if information is provided detector level effects
+ if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
+ if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
Bool_t selectEvent = false;
Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
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 take total momentum here
- fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
- jInp.set_user_index(i);
- inputParticlesRec.push_back(jInp);
+
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
+ //Add particles to fastjet in case we are not running toy model
+ fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->P());
+ jInp.set_user_index(i);
+ inputParticlesRec.push_back(jInp);
+ }
+ else if(fUseDiceEfficiency) {
+
+ // Dice to decide if particle is kept or not - toy model for efficiency
+ //
+ Double_t rnd = fRandom->Uniform(1.);
+ Double_t pT = vp->Pt();
+ Double_t eff[3] = {0.};
+ Double_t pTtmp = pT;
+ if(pT>10.) pTtmp = 10.;
+ Double_t eff1 = fhEffH1->GetBinContent(fhEffH1->FindBin(pTtmp));
+ Double_t eff2 = fhEffH2->GetBinContent(fhEffH2->FindBin(pTtmp));
+ Double_t eff3 = fhEffH3->GetBinContent(fhEffH3->FindBin(pTtmp));
+ Int_t cat[3] = {0};
+ //Sort efficiencies from large to small
+ if(eff1>eff2 && eff1>eff3) {
+ eff[0] = eff1;
+ cat[0] = 1;
+ if(eff2>eff3) {
+ eff[1] = eff2;
+ eff[2] = eff3;
+ cat[1] = 2;
+ cat[2] = 3;
+ } else {
+ eff[1] = eff3;
+ eff[2] = eff2;
+ cat[1] = 3;
+ cat[2] = 2;
+ }
+ }
+ else if(eff2>eff1 && eff2>eff3) {
+ eff[0] = eff2;
+ cat[0] = 2;
+ if(eff1>eff3) {
+ eff[1] = eff1;
+ eff[2] = eff3;
+ cat[1] = 1;
+ cat[2] = 3;
+ } else {
+ eff[1] = eff3;
+ eff[2] = eff1;
+ cat[1] = 3;
+ cat[2] = 1;
+ }
+ }
+ else if(eff3>eff1 && eff3>eff2) {
+ eff[0] = eff3;
+ cat[0] = 3;
+ if(eff1>eff2) {
+ eff[1] = eff1;
+ eff[2] = eff2;
+ cat[1] = 1;
+ cat[2] = 2;
+ } else {
+ eff[1] = eff2;
+ eff[2] = eff1;
+ cat[1] = 2;
+ cat[2] = 1;
+ }
+ }
+
+ Double_t sumEff = eff[0]+eff[1]+eff[2];
+ fp1Efficiency->Fill(vp->Pt(),sumEff);
+ if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
+
+ if(fUseTrPtResolutionSmearing) {
+ //Smear momentum of generated particle
+ Double_t smear = 1.;
+ //Select hybrid track category
+ if(rnd<=eff[2])
+ smear = GetMomentumSmearing(cat[2],pT);
+ else if(rnd<=(eff[2]+eff[1]))
+ smear = GetMomentumSmearing(cat[1],pT);
+ else
+ smear = GetMomentumSmearing(cat[0],pT);
+
+ fp1PtResolution->Fill(vp->Pt(),smear);
+
+ Double_t sigma = vp->Pt()*smear;
+ Double_t pTrec = fRandom->Gaus(vp->Pt(),sigma);
+
+ Double_t phi = vp->Phi();
+ Double_t theta = 2.*TMath::ATan(TMath::Exp(-1.*(vp->Eta())));
+ Double_t pX = pTrec * TMath::Cos(phi);
+ Double_t pY = pTrec * TMath::Sin(phi);
+ Double_t pZ = pTrec/TMath::Tan(theta);
+ Double_t p=TMath::Sqrt(pTrec*pTrec+pZ*pZ);
+
+ fh2PtGenPtSmeared->Fill(vp->Pt(),pTrec);
+
+ fastjet::PseudoJet jInp(pX,pY,pZ,p);
+ jInp.set_user_index(i);
+ inputParticlesRec.push_back(jInp);
+
+ }
+ else {
+ 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
if(fTCAJetsOut){
if(i == 0){
fRef->Delete(); // make sure to delete before placement new...
- new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
+ new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
+ }
}
- fRef->Add(vp);
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
}
}// recparticles
// 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);
Float_t pTback = 0;
if(externalBackground){
// carefull has to be filled in a task before
- // todo, ReArrange to the botom
- pTback = externalBackground->GetBackground(1)*leadingJet.EffectiveAreaCharged();
+ // todo, ReArrange to the botom
+ pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged();
}
pt = leadingJet.Pt() - pTback;
// correlation of leading jet with tracks
fh2NConstPt->Fill(tmpPt,constituents.size());
// loop over constiutents and fill spectrum
+ AliVParticle *partLead = 0;
+ Float_t ptLead = -1;
+
for(unsigned int ic = 0; ic < constituents.size();ic++){
AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
+ if(!part) continue;
fh1PtJetConstRec->Fill(part->Pt());
if(aodOutJet){
- aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
+ if(part->Pt()>fMaxTrackPtInJet){
+ aodOutJet->SetTrigger(AliAODJet::kHighTrackPtTriggered);
+ }
+ }
+ if(part->Pt()>ptLead){
+ ptLead = part->Pt();
+ partLead = part;
}
if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
}
-
+
+ AliAODTrack *aodT = 0;
+ if(partLead){
+ if(aodOutJet){
+ //set pT of leading constituent of jet
+ aodOutJet->SetPtLeading(partLead->Pt());
+ aodT = dynamic_cast<AliAODTrack*>(partLead);
+ if(aodT){
+ if(aodT->TestFilterBit(fFilterMaskBestPt)){
+ aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
+ }
+ }
+ }
+ }
+
// correlation
Float_t tmpPhi = tmpRec.Phi();
Float_t tmpEta = tmpRec.Eta();
}
if(skip)continue;
tmpRecC.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below
+ tmpRecC.SetPtLeading(-1.);
if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC);
if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC);
}// loop over random cones creation
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);
+ if(vp->Pt() > jC->GetPtLeading()) jC->SetPtLeading(vp->Pt());
}
}
}// 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();
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);
+ if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt());
}
}
}
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
+ // rescale the momentum vectors for the random cones
AliAODJet *rC = (AliAODJet*)fTCARandomConesOut->At(ir);
if(rC){
}// if(fNRandomCones
//background estimates:all bckg jets(0) & wo the 2 hardest(1)
-
-
-
-
if(fAODJetBackgroundOut){
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.;
}
}
+
+ if(fStoreRhoLeadingTrackCorr) {
+ vector<fastjet::PseudoJet> jets3=sortedJets;
+ if(jets3.size()>2) jets3.erase(jets3.begin(),jets3.begin()+2);
+
+ Double_t bkg2=0;
+ Double_t sigma2=0.;
+ Double_t meanarea2=0.;
+
+ clustSeq.get_median_rho_and_sigma(jets3, range, false, bkg2, sigma2, meanarea2, true);
+ fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2);
+
+ fh2CentvsRho->Fill(cent,bkg2);
+ fh2CentvsSigma->Fill(cent,sigma2);
+ fh2MultvsRho->Fill(nCh,bkg2);
+ fh2MultvsSigma->Fill(nCh,sigma2);
+
+
+ //Calculate rho with 4-vector area method for different quadrants with respect to the leading track in the event
+ //exclude 2 leading jets from event
+ //Quadrant 1: |phi_leadingTrack - phi_bkgCluster| < pi/2./2. - R (Near side to leading track)
+ //Quadrant 2: pi/2 - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < pi/2 + (pi/2./2. - R) (Orthogonal to leading track)
+ //Quadrant 3: pi - (pi/2/2 - R) < |phi_leadingTrack - phi_bkgCluster| < pi + (pi/2./2. - R) (Away side to leading track)
+ //Quadrant 4: 3/2*pi - (pi/2./2. - R) < |phi_leadingTrack - phi_bkgCluster| < 3/2*pi + (pi/2./2. - R) (Orthogonal to leading track)
+
+ //Assuming R=0.2 for background clusters
+
+ Double_t bkg2Q[4] = {0.};
+ Double_t sigma2Q[4] = {0.};
+ Double_t meanarea2Q[4] = {0.};
+
+ //Get phi of leading track in event
+ AliVParticle *leading = (AliVParticle*)recParticles.At(0);
+ Float_t phiLeadingTrack = leading->Phi();
+ Float_t phiStep = (TMath::Pi()/2./2. - 0.2);
+
+ //Quadrant1 - near side
+ phiMin = phiLeadingTrack - phiStep;
+ phiMax = phiLeadingTrack + phiStep;
+ fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets3, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
+
+ //Quadrant2 - orthogonal
+ Double_t phiQ2 = phiLeadingTrack + TMath::Pi()/2.;
+ if(phiQ2>TMath::TwoPi()) phiQ2 = phiQ2 - TMath::TwoPi();
+ phiMin = phiQ2 - phiStep;
+ phiMax = phiQ2 + phiStep;
+ fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets3, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
+
+ //Quadrant3 - away side
+ Double_t phiQ3 = phiLeadingTrack + TMath::Pi();
+ if(phiQ3>TMath::TwoPi()) phiQ3 = phiQ3 - TMath::TwoPi();
+ phiMin = phiQ3 - phiStep;
+ phiMax = phiQ3 + phiStep;
+ fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets3, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
+
+ //Quadrant4 - othogonal
+ Double_t phiQ4 = phiLeadingTrack + 3./2.*TMath::Pi();
+ if(phiQ4>TMath::TwoPi()) phiQ4 = phiQ4 - TMath::TwoPi();
+ phiMin = phiQ4 - phiStep;
+ phiMax = phiQ4 + phiStep;
+ fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets3, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
+
+ fh3CentvsRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0],leading->Pt());
+ fh3CentvsRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1],leading->Pt());
+ fh3CentvsRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2],leading->Pt());
+ fh3CentvsRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3],leading->Pt());
+
+ fh3CentvsSigmaLeadingTrackPtQ1->Fill(cent,sigma2Q[0],leading->Pt());
+ fh3CentvsSigmaLeadingTrackPtQ2->Fill(cent,sigma2Q[1],leading->Pt());
+ fh3CentvsSigmaLeadingTrackPtQ3->Fill(cent,sigma2Q[2],leading->Pt());
+ fh3CentvsSigmaLeadingTrackPtQ4->Fill(cent,sigma2Q[3],leading->Pt());
+
+ fh3MultvsRhoLeadingTrackPtQ1->Fill(nCh,bkg2Q[0],leading->Pt());
+ fh3MultvsRhoLeadingTrackPtQ2->Fill(nCh,bkg2Q[1],leading->Pt());
+ fh3MultvsRhoLeadingTrackPtQ3->Fill(nCh,bkg2Q[2],leading->Pt());
+ fh3MultvsRhoLeadingTrackPtQ4->Fill(nCh,bkg2Q[3],leading->Pt());
+
+ fh3MultvsSigmaLeadingTrackPtQ1->Fill(nCh,sigma2Q[0],leading->Pt());
+ fh3MultvsSigmaLeadingTrackPtQ2->Fill(nCh,sigma2Q[1],leading->Pt());
+ fh3MultvsSigmaLeadingTrackPtQ3->Fill(nCh,sigma2Q[2],leading->Pt());
+ fh3MultvsSigmaLeadingTrackPtQ4->Fill(nCh,sigma2Q[3],leading->Pt());
+
+ }
void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/)
{
-// Terminate analysis
-//
+ //
+ // Terminate analysis
+ //
if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n");
+
+ if(fMomResH1Fit) delete fMomResH1Fit;
+ if(fMomResH2Fit) delete fMomResH2Fit;
+ if(fMomResH3Fit) delete fMomResH3Fit;
+
}
Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){
+ //
+ // get list of tracks/particles for different types
+ //
+
if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
Int_t iCount = 0;
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->IsHybridITSTPC();
- else if(fFilterType == 2)bGood = tr->IsHybridTPC();
+ 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(!trackAOD)continue;
Bool_t bGood = false;
if(fFilterType == 0)bGood = true;
- else if(fFilterType == 1)bGood = trackAOD->IsHybridITSTPC();
- else if(fFilterType == 2)bGood = trackAOD->IsHybridTPC();
+ 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;
return iCount;
}
+void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() {
+
+ TFile *f = TFile::Open(fPathTrPtResolution.Data());
+
+ if(!f)return;
+
+ TProfile *fProfPtPtSigma1PtGlobSt = (TProfile*)f->Get("fProfPtPtSigma1PtGlobSt");
+ TProfile *fProfPtPtSigma1PtGlobCnoITS = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoITS");
+ TProfile *fProfPtPtSigma1PtGlobCnoSPD = (TProfile*)f->Get("fProfPtPtSigma1PtGlobCnoSPD");
+
+ SetSmearResolution(kTRUE);
+ SetMomentumResolutionHybrid(fProfPtPtSigma1PtGlobSt,fProfPtPtSigma1PtGlobCnoITS,fProfPtPtSigma1PtGlobCnoSPD);
+
+
+}
+
+void AliAnalysisTaskJetCluster::LoadTrEfficiencyRootFileFromOADB() {
+
+ TFile *f = TFile::Open(fPathTrEfficiency.Data());
+ if(!f)return;
+
+ TH1D *hEffPosGlobSt = (TH1D*)f->Get("hEffPosGlobSt");
+ TH1D *hEffPosGlobCnoITS = (TH1D*)f->Get("hEffPosGlobCnoITS");
+ TH1D *hEffPosGlobCnoSPD = (TH1D*)f->Get("hEffPosGlobCnoSPD");
+
+ SetDiceEfficiency(kTRUE);
+
+ if(fChangeEfficiencyFraction>0.) {
+
+ TH1D *hEffPosGlobStMin = (TH1D*)hEffPosGlobSt->Clone("hEffPosGlobStMin");
+
+ for(int bin=1; bin<=hEffPosGlobSt->GetNbinsX(); bin++) {
+ Double_t content = hEffPosGlobSt->GetBinContent(bin);
+ hEffPosGlobStMin->SetBinContent(bin,content-fChangeEfficiencyFraction);
+ }
+
+ SetEfficiencyHybrid(hEffPosGlobStMin,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
+
+ }
+ else
+ SetEfficiencyHybrid(hEffPosGlobSt,hEffPosGlobCnoITS,hEffPosGlobCnoSPD);
+
+}
+
+void AliAnalysisTaskJetCluster::SetMomentumResolutionHybrid(TProfile *p1, TProfile *p2, TProfile *p3) {
+
+ //
+ // set mom res profiles
+ //
+
+ if(fMomResH1) delete fMomResH1;
+ if(fMomResH2) delete fMomResH2;
+ if(fMomResH3) delete fMomResH3;
+
+ fMomResH1 = new TProfile(*p1);//(TProfile*)p1->Clone("fMomResH1");
+ fMomResH2 = new TProfile(*p2);//(TProfile*)p2->Clone("fMomResH2");
+ fMomResH3 = new TProfile(*p3);//(TProfile*)p3->Clone("fMomResH3");
+
+}
+
+void AliAnalysisTaskJetCluster:: SetEfficiencyHybrid(TH1 *h1, TH1 *h2, TH1 *h3) {
+ //
+ // set tracking efficiency histos
+ //
+
+ fhEffH1 = (TH1*)h1->Clone("fhEffH1");
+ fhEffH2 = (TH1*)h2->Clone("fhEffH2");
+ fhEffH3 = (TH1*)h3->Clone("fhEffH3");
+}
+
+Double_t AliAnalysisTaskJetCluster::GetMomentumSmearing(Int_t cat, Double_t pt) {
+
+ //
+ // Get smearing on generated momentum
+ //
+
+ //printf("GetMomentumSmearing for cat %d and pt = %f \n",cat,pt);
+
+ TProfile *fMomRes = 0x0;
+ if(cat==1) fMomRes = (TProfile*)fMomResH1->Clone("fMomRes");
+ if(cat==2) fMomRes = (TProfile*)fMomResH2->Clone("fMomRes");
+ if(cat==3) fMomRes = (TProfile*)fMomResH3->Clone("fMomRes");
+
+ if(!fMomRes) {
+ return 0.;
+ }
+
+
+ Double_t smear = 0.;
+
+ if(pt>20.) {
+ if(cat==1 && fMomResH1Fit) smear = fMomResH1Fit->Eval(pt);
+ if(cat==2 && fMomResH2Fit) smear = fMomResH2Fit->Eval(pt);
+ if(cat==3 && fMomResH3Fit) smear = fMomResH3Fit->Eval(pt);
+ }
+ else {
+
+ Int_t bin = fMomRes->FindBin(pt);
+
+ smear = fRandom->Gaus(fMomRes->GetBinContent(bin),fMomRes->GetBinError(bin));
+
+ }
+
+ if(fMomRes) delete fMomRes;
+
+ return smear;
+}
+
+void AliAnalysisTaskJetCluster::FitMomentumResolution() {
+ //
+ // Fit linear function on momentum resolution at high pT
+ //
+
+ if(!fMomResH1Fit && fMomResH1) {
+ fMomResH1Fit = new TF1("fMomResH1Fit","[0]+[1]*x",0.,200.);
+ fMomResH1->Fit(fMomResH1Fit,"LL V0","",5.,30.);
+ fMomResH1Fit ->SetRange(5.,100.);
+ }
+
+ if(!fMomResH2Fit && fMomResH2) {
+ fMomResH2Fit = new TF1("fMomResH2Fit","[0]+[1]*x",0.,200.);
+ fMomResH2->Fit(fMomResH2Fit,"LL V0","",5.,30.);
+ fMomResH2Fit ->SetRange(5.,100.);
+ }
+
+ if(!fMomResH3Fit && fMomResH3) {
+ fMomResH3Fit = new TF1("fMomResH3Fit","[0]+[1]*x",0.,200.);
+ fMomResH3->Fit(fMomResH3Fit,"LL V0","",5.,30.);
+ fMomResH3Fit ->SetRange(5.,100.);
+ }
+
+}
+
/*
Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector<fastjet::PseudoJet> &inputParticles){
for(int i = 0; i < particles.GetEntries(); i++){
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff