// get info on how fastjet was configured
#include "fastjet/config.h"
+using std::vector;
ClassImp(AliAnalysisTaskJetCluster)
fhEffH1(0x0),
fhEffH2(0x0),
fhEffH3(0x0),
- fUseTrMomentumSmearing(kFALSE),
+ 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),
+ fh3CentvsRhoLeadingTrackCorr(0x0),
+ fh3CentvsSigmaLeadingTrackCorr(0x0),
+ fh3MultvsRhoLeadingTrackCorr(0x0),
+ fh3MultvsSigmaLeadingTrackCorr(0x0),
fh2PtGenPtSmeared(0x0),
fp1Efficiency(0x0),
fp1PtResolution(0x0),
fhEffH1(0x0),
fhEffH2(0x0),
fhEffH3(0x0),
- fUseTrMomentumSmearing(kFALSE),
+ 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),
+ fh3CentvsRhoLeadingTrackCorr(0x0),
+ fh3CentvsSigmaLeadingTrackCorr(0x0),
+ fh3MultvsRhoLeadingTrackCorr(0x0),
+ fh3MultvsSigmaLeadingTrackCorr(0x0),
fh2PtGenPtSmeared(0x0),
fp1Efficiency(0x0),
fp1PtResolution(0x0),
if(fJetTypes&kJetRan){
fTCAJetsOutRan = new TClonesArray("AliAODJet", 0);
fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random"));
- if(fUseDiceEfficiency || fUseTrMomentumSmearing)
- fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%d",fNonStdBranch.Data(),"random",fUseTrMomentumSmearing,fUseDiceEfficiency));
+ 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(!AODEvent()->FindListObject(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()))){
fAODJetBackgroundOut = new AliAODJetEventBackground();
fAODJetBackgroundOut->SetName(Form("%s_%s",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data()));
- if(fUseDiceEfficiency || fUseTrMomentumSmearing)
- fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency));
+ 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 || fUseTrMomentumSmearing)
- cName = Form("%sDetector%d%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrMomentumSmearing,fUseDiceEfficiency,fNSkipLeadingCone);
+ if(fUseDiceEfficiency || fUseTrPtResolutionSmearing)
+ cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone);
if(fNRandomCones>0){
if(fJetTypes&kRC){
}
}
- // FitMomentumResolution();
-
if(!fHistList)fHistList = new TList();
fHistList->SetOwner();
fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
+ fh3CentvsRhoLeadingTrackCorr = new TH3F("fh3CentvsRhoLeadingTrackCorr","centrality vs background density; centrality; #rho; Quadrant", 100,0.,100.,600,0.,300.,4,0.,4.);
+ fh3CentvsSigmaLeadingTrackCorr = new TH3F("fh3CentvsSigmaLeadingTrackCorr","centrality vs backgroun sigma; centrality; #sigma(#rho); Quadrant",100,0.,100.,500,0.,50.,4,0.,4.);
+ fHistList->Add(fh3CentvsRhoLeadingTrackCorr);
+ fHistList->Add(fh3CentvsSigmaLeadingTrackCorr);
+
+ fh3MultvsRhoLeadingTrackCorr = new TH3F("fh3MultvsRhoLeadingTrackCorr","mult vs background density; N_{tracks}; #rho; Quadrant", 100,0.,5000.,600,0.,300.,4,0.,4.);
+ fh3MultvsSigmaLeadingTrackCorr = new TH3F("fh3MultvsSigmaLeadingTrackCorr","mult vs background sigma; N_{tracks}; #sigma(#rho); Quadrant",100,0.,5000.,500,0.,50.,4,0.,4.);
+ fHistList->Add(fh3MultvsRhoLeadingTrackCorr);
+ fHistList->Add(fh3MultvsSigmaLeadingTrackCorr);
+
//Detector level effects histos
fh2PtGenPtSmeared = new TH2F("fh2PtGenPtSmeared","fh2PtGenPtSmeared",nBinPt,binLimitsPt,nBinPt,binLimitsPt);
TH1::AddDirectory(oldStatus);
}
-void AliAnalysisTaskJetCluster::Init()
+void AliAnalysisTaskJetCluster::LocalInit()
{
//
// Initialization
if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n");
+ if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB();
+ if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB();
+
+
FitMomentumResolution();
}
}
//Check if information is provided detector level effects
- if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrMomentumSmearing = kFALSE;
+ if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE;
if(!fhEffH1 || !fhEffH2 || !fhEffH3) fUseDiceEfficiency = kFALSE;
Bool_t selectEvent = false;
// Carefull energy is not well determined in real data, should not matter for p_T scheme?
// we take total momentum here
- if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
+ 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);
Double_t sumEff = eff[0]+eff[1]+eff[2];
fp1Efficiency->Fill(vp->Pt(),sumEff);
- if(rnd>sumEff) continue;
+ if(rnd>sumEff && pT > fDiceEfficiencyMinPt) continue;
- if(fUseTrMomentumSmearing) {
+ if(fUseTrPtResolutionSmearing) {
//Smear momentum of generated particle
Double_t smear = 1.;
//Select hybrid track category
if(fTCAJetsOut){
if(i == 0){
fRef->Delete(); // make sure to delete before placement new...
- if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) {
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) {
new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp)); //TRefArray does not work with toy model ...
}
}
- if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
+ if((!fUseTrPtResolutionSmearing) && (!fUseDiceEfficiency)) fRef->Add(vp); //TRefArray does not work with toy model ...
}
}// recparticles
if(!part) continue;
fh1PtJetConstRec->Fill(part->Pt());
if(aodOutJet){
- if((!fUseTrMomentumSmearing) && (!fUseDiceEfficiency)) 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){
- aodT = dynamic_cast<AliAODTrack*>(partLead);
- if(aodT){
- if(aodT->TestFilterBit(fFilterMaskBestPt)){
- aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest);
+ 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
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();
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){
// fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone));
// fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone));
+
+ //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();
+
+ //Quadrant1 - near side
+ phiMin = phiLeadingTrack - (TMath::Pi()/2./2. - 0.2);
+ phiMax = phiLeadingTrack + (TMath::Pi()/2./2. - 0.2);
+ fastjet::RangeDefinition rangeQ1(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets2, rangeQ1, false, bkg2Q[0], sigma2Q[0], meanarea2Q[0], true);
+
+ //Quadrant2 - orthogonal
+ phiMin = phiLeadingTrack + TMath::Pi()/2. - (TMath::Pi()/2./2. - 0.2);
+ phiMax = phiLeadingTrack + TMath::Pi()/2. + (TMath::Pi()/2./2. - 0.2);
+ fastjet::RangeDefinition rangeQ2(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets2, rangeQ2, false, bkg2Q[1], sigma2Q[1], meanarea2Q[1], true);
+
+ //Quadrant3 - away side
+ phiMin = phiLeadingTrack + TMath::Pi() - (TMath::Pi()/2./2. - 0.2);
+ phiMax = phiLeadingTrack + TMath::Pi() + (TMath::Pi()/2./2. - 0.2);
+ fastjet::RangeDefinition rangeQ3(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets2, rangeQ3, false, bkg2Q[2], sigma2Q[2], meanarea2Q[2], true);
+
+ //Quadrant4 - othogonal
+ phiMin = phiLeadingTrack + TMath::Pi()*3./2. - (TMath::Pi()/2./2. - 0.2);
+ phiMax = phiLeadingTrack + TMath::Pi()*3./2. + (TMath::Pi()/2./2. - 0.2);
+ fastjet::RangeDefinition rangeQ4(rapMin,rapMax, phiMin, phiMax);
+ clustSeq.get_median_rho_and_sigma(jets2, rangeQ4, false, bkg2Q[3], sigma2Q[3], meanarea2Q[3], true);
+
+ fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[0],0.5);
+ fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[1],1.5);
+ fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[2],2.5);
+ fh3CentvsRhoLeadingTrackCorr->Fill(cent,bkg2Q[3],3.5);
+
+ fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[0],0.5);
+ fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[1],1.5);
+ fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[2],2.5);
+ fh3CentvsSigmaLeadingTrackCorr->Fill(cent,sigma2Q[3],3.5);
+
+ fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[0],0.5);
+ fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[1],1.5);
+ fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[2],2.5);
+ fh3MultvsRhoLeadingTrackCorr->Fill(nCh,bkg2Q[3],3.5);
+
+ fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[0],0.5);
+ fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[1],1.5);
+ fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[2],2.5);
+ fh3MultvsSigmaLeadingTrackCorr->Fill(nCh,sigma2Q[3],3.5);
+
}
}
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
//
- fMomResH1 = (TProfile*)p1->Clone("fMomResH1");
- fMomResH2 = (TProfile*)p2->Clone("fMomResH2");
- fMomResH3 = (TProfile*)p3->Clone("fMomResH3");
+ 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) {
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff