X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=JETAN%2FAliAnalysisTaskJetCluster.cxx;h=8173f51041f1c3a62c4f2234669ba11073332811;hb=700a5a265e0fc00f9fec37f5d0fb5fed00d81c6b;hp=f4422b9ba2654c990bdf7b7e156e99f4eb810d5d;hpb=0341d0d825a887ac2507fdef49b46df5c9e51606;p=u%2Fmrichter%2FAliRoot.git diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index f4422b9ba26..8173f51041f 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -32,10 +32,12 @@ #include #include #include +#include #include #include #include #include "TDatabasePDG.h" +#include #include "AliAnalysisTaskJetCluster.h" #include "AliAnalysisManager.h" @@ -45,6 +47,7 @@ #include "AliESDEvent.h" #include "AliAODEvent.h" #include "AliAODHandler.h" +#include "AliAODExtension.h" #include "AliAODTrack.h" #include "AliAODJet.h" #include "AliAODMCParticle.h" @@ -64,36 +67,91 @@ // 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; } -AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), +AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): + AliAnalysisTaskSE(), fAOD(0x0), fAODExtension(0x0), fRef(new TRefArray), fUseAODTrackInput(kFALSE), fUseAODMCInput(kFALSE), - fUseGlobalSelection(kFALSE), fUseBackgroundCalc(kFALSE), + fEventSelection(kFALSE), + fRequireVZEROAC(kFALSE), + fRequireTZEROvtx(kFALSE), + fUseHFcuts(kFALSE), fFilterMask(0), + fFilterMaskBestPt(0), + fFilterType(0), + fJetTypes(kJet), fTrackTypeRec(kTrackUndef), fTrackTypeGen(kTrackUndef), fNSkipLeadingRan(0), - fNRandomCones(5), + fNSkipLeadingCone(0), + fNRandomCones(0), fAvgTrials(1), - fExternalWeight(1), + fExternalWeight(1), + fTrackEtaWindow(0.9), + fRequireITSRefit(0), + fApplySharedClusterCut(0), fRecEtaWindow(0.5), fTrackPtCut(0.), - fJetOutputMinPt(1), + 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(-1.), + fUseTrPtResolutionFromOADB(kFALSE), + fUseTrEfficiencyFromOADB(kFALSE), + fPathTrPtResolution(""), + fPathTrEfficiency(""), + fChangeEfficiencyFraction(0.), + fEfficiencyFixed(1.), fRparam(1.0), fAlgorithm(fastjet::kt_algorithm), fStrategy(fastjet::Best), @@ -102,6 +160,11 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fGhostArea(0.01), fActiveAreaRepeats(1), fGhostEtamax(1.5), + fTCAJetsOut(0x0), + fTCAJetsOutRan(0x0), + fTCARandomConesOut(0x0), + fTCARandomConesOutRan(0x0), + fAODJetBackgroundOut(0x0), fRandom(0), fh1Xsec(0x0), fh1Trials(0x0), @@ -124,8 +187,12 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fh1PtJetsRecInRan(0x0), fh1PtTracksGenIn(0x0), fh1Nch(0x0), + fh1CentralityPhySel(0x0), fh1Centrality(0x0), fh1CentralitySelect(0x0), + fh1ZPhySel(0x0), + fh1Z(0x0), + fh1ZSelect(0x0), fh2NRecJetsPt(0x0), fh2NRecTracksPt(0x0), fh2NConstPt(0x0), @@ -153,12 +220,49 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(), fh2PtNchNRan(0x0), fh2TracksLeadingJetPhiPtRan(0x0), fh2TracksLeadingJetPhiPtWRan(0x0), + fh3CentvsRhoLeadingTrackPt(0x0), + fh3CentvsSigmaLeadingTrackPt(0x0), + fh3MultvsRhoLeadingTrackPt(0x0), + fh3MultvsSigmaLeadingTrackPt(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), + fh3CentvsDeltaRhoLeadingTrackPtQ1(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ2(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ3(0x0), + fh3CentvsDeltaRhoLeadingTrackPtQ4(0x0), + fh2PtGenPtSmeared(0x0), + fp1Efficiency(0x0), + fp1PtResolution(0x0), fHistList(0x0) { + // + // Constructor + // + for(int i = 0;i<3;i++){ fh1BiARandomCones[i] = 0; fh1BiARandomConesRan[i] = 0; } + for(int i = 0;i cleared in the UserExec.... // here we can also have the case that the brnaches are written to a separate file + + if(fJetTypes&kJet){ + fTCAJetsOut = new TClonesArray("AliAODJet", 0); + fTCAJetsOut->SetName(fNonStdBranch.Data()); + AddAODBranch("TClonesArray",&fTCAJetsOut,fNonStdFile.Data()); + } - TClonesArray *tca = new TClonesArray("AliAODJet", 0); - tca->SetName(fNonStdBranch.Data()); - AddAODBranch("TClonesArray",&tca,fNonStdFile.Data()); - - - TClonesArray *tcaran = new TClonesArray("AliAODJet", 0); - tcaran->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); - AddAODBranch("TClonesArray",&tcaran,fNonStdFile.Data()); - + if(fJetTypes&kJetRan){ + fTCAJetsOutRan = new TClonesArray("AliAODJet", 0); + fTCAJetsOutRan->SetName(Form("%s_%s",fNonStdBranch.Data(),"random")); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + fTCAJetsOutRan->SetName(Form("%s_%sDetector%d%dEffFixed%d",fNonStdBranch.Data(),"random",fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.))); + else + 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()))){ - 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())); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + fAODJetBackgroundOut->SetName(Form("%s_%sDetector%d%dEffFixed%d",AliAODJetEventBackground::StdBranchName(),fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.))); + else + 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("%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 for the random cones with the same R + TString cName = Form("%sRandomConeSkip%02d",fNonStdBranch.Data(),fNSkipLeadingCone); + if(fUseDiceEfficiency || fUseTrPtResolutionSmearing) { + if( fEfficiencyFixed < 1.) + cName = Form("%sDetector%d%dEffFixed%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fEfficiencyFixed*100.),fNSkipLeadingCone); + else + cName = Form("%sDetector%d%dFr%d_RandomConeSkip%02d",fNonStdBranch.Data(),fUseTrPtResolutionSmearing,fUseDiceEfficiency,(int)(fChangeEfficiencyFraction*100.),fNSkipLeadingCone); + } + if(fNRandomCones>0){ + 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(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); @@ -354,14 +550,14 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() // // 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; } } @@ -389,7 +585,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() } } - const int nChMax = 4000; + const int nChMax = 5000; fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); @@ -417,13 +613,18 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() 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); 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); @@ -454,43 +655,124 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2TrackEtaPt = new TH2F("fh2TrackEtaPt","pt vs eta all jets;#eta;p_{T}", - nBinEta,binLimitsEta,nBinPt,binLimitsPt); + nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2LeadingTrackEtaPt = new TH2F("fh2LeadingTrackEtaPt","pT vs eta leading jets;#eta;p_{T}", - nBinEta,binLimitsEta,nBinPt,binLimitsPt); + nBinEta,binLimitsEta,nBinPt,binLimitsPt); fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta", - nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); + nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta", nBinPhi,binLimitsPhi,nBinEta,binLimitsEta); fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPt = new TH2F("fh2TracksLeadingJetPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtRan = new TH2F("fh2TracksLeadingJetPhiPtRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2JetsLeadingPhiPtW = new TH2F("fh2JetsLeadingPhiPtW","leading p_T vs delta phi p_T weigted to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingPhiPtW = new TH2F("fh2TracksLeadingPhiPtW","leading p_T vs delta phi to leading jet (p_T weighted);#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtW = new TH2F("fh2TracksLeadingJetPhiPtW","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); fh2TracksLeadingJetPhiPtWRan = new TH2F("fh2TracksLeadingJetPhiPtWRan","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)", - nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); + + if(fStoreRhoLeadingTrackCorr) { + fh3CentvsRhoLeadingTrackPt = new TH3F("fh3CentvsRhoLeadingTrackPt","centrality vs background density full event; centrality; #rho", 50,0.,100.,500,0.,250.,100,0.,100.); + fh3CentvsSigmaLeadingTrackPt = new TH3F("fh3CentvsSigmaLeadingTrackPt","centrality vs sigma full event; centrality; #sigma(#rho)", 50,0.,100.,50,0.,50.,100,0.,100.); + fh3MultvsRhoLeadingTrackPt = new TH3F("fh3MultvsRhoLeadingTrackPt","multiplicity vs background density full event; multiplicity; #rho", 100,0.,5000.,500,0.,250.,100,0.,100.); + fh3MultvsSigmaLeadingTrackPt = new TH3F("fh3MultvsSigmaLeadingTrackPt","multiplicity vs sigma full event; multiplicity; #sigma(#rho)", 100,0.,5000.,50,0.,50.,100,0.,100.); + + + 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.); + + + fh3CentvsDeltaRhoLeadingTrackPtQ1 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ1","centrality vs background density Q1; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ2 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ2","centrality vs background density Q2; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ3 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ3","centrality vs background density Q3; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + fh3CentvsDeltaRhoLeadingTrackPtQ4 = new TH3F("fh3CentvsDeltaRhoLeadingTrackPtQ4","centrality vs background density Q4; centrality; #rho_{quadrant}-#rho_{full event}; leading p_{t}^{track}", 50,0.,100.,200,-10.,10.,100,0.,100.); + + fHistList->Add(fh3CentvsRhoLeadingTrackPt); + fHistList->Add(fh3CentvsSigmaLeadingTrackPt); + fHistList->Add(fh3MultvsRhoLeadingTrackPt); + fHistList->Add(fh3MultvsSigmaLeadingTrackPt); + + 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); + + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ1); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ2); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ3); + fHistList->Add(fh3CentvsDeltaRhoLeadingTrackPtQ4); + } + + //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); - if(fUseBackgroundCalc){ + fHistList->Add(fh2PtGenPtSmeared); + fHistList->Add(fp1Efficiency); + fHistList->Add(fp1PtResolution); + + if(fNRandomCones>0&&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(fh1Xsec); @@ -513,12 +795,23 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh1Nch); fHistList->Add(fh1Centrality); fHistList->Add(fh1CentralitySelect); - if(fUseBackgroundCalc){ + fHistList->Add(fh1CentralityPhySel); + fHistList->Add(fh1Z); + fHistList->Add(fh1ZSelect); + fHistList->Add(fh1ZPhySel); + if(fNRandomCones>0&&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); @@ -546,7 +839,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() fHistList->Add(fh2NConstLeadingPtRan); fHistList->Add(fh2TracksLeadingJetPhiPtRan); fHistList->Add(fh2TracksLeadingJetPhiPtWRan); - } + } // =========== Switch on Sumw2 for all histos =========== for (Int_t i=0; iGetEntries(); ++i) { @@ -561,7 +854,7 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects() TH1::AddDirectory(oldStatus); } -void AliAnalysisTaskJetCluster::Init() +void AliAnalysisTaskJetCluster::LocalInit() { // // Initialization @@ -569,55 +862,30 @@ void AliAnalysisTaskJetCluster::Init() if (fDebug > 1) printf("AnalysisTaskJetCluster::Init() \n"); -} + if(fUseTrPtResolutionFromOADB) LoadTrPtResolutionRootFileFromOADB(); + if(fUseTrEfficiencyFromOADB) LoadTrEfficiencyRootFileFromOADB(); -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; - } + FitMomentumResolution(); +} +void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) +{ // 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(); - - 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(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());; } // // Execute analysis for current event @@ -629,7 +897,7 @@ 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... @@ -642,29 +910,87 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fESD = dynamic_cast (InputEvent()); } } - + + //Check if information is provided detector level effects + if(!fMomResH1 || !fMomResH2 || !fMomResH3) fUseTrPtResolutionSmearing = kFALSE; + if( fEfficiencyFixed < 1. ) { + if (!fUseDiceEfficiency) + fUseDiceEfficiency = 1; // 1 is the default; 2 can be set by user + } + else { + if(!fhEffH1 || !fhEffH2 || !fhEffH3 ) fUseDiceEfficiency = kFALSE; + } + 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(); + 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)<20.&&r2<1.){ // apply vertex cut later on - if(physicsSelection)selectEvent = true; + if(TMath::Abs(zVtx)0){ + if(centfCentCutUp){ + selectEvent = false; } + } + }else{ + selectEvent = true; } } + + + Bool_t T0 = false; + Bool_t V0 = false; + const AliAODVZERO *vzero = fAOD->GetVZEROData(); + if(vzero){ + if((vzero->GetTriggerChargeA()>0)&&(vzero->GetTriggerChargeC()>0)){ + V0 = true; + } + } + + const AliAODTZERO *tzero = fAOD->GetTZEROData(); + if(tzero){ + if(TMath::Abs(tzero->GetT0VertexRaw())<100){ + T0 = true; + } + } + + if(fRequireVZEROAC&&fRequireTZEROvtx)selectEvent = selectEvent&&V0&&T0; + else if(fRequireTZEROvtx)selectEvent = selectEvent&&T0; + else if(fRequireVZEROAC)selectEvent = selectEvent&&V0; + + if(!selectEvent){ PostData(1, fHistList); return; } - + fh1Centrality->Fill(cent); + fh1Z->Fill(zVtx); fh1Trials->Fill("#sum{ntrials}",1); @@ -681,6 +1007,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec); Float_t nCh = recParticles.GetEntries(); + Float_t nGen=genParticles.GetEntries(); fh1Nch->Fill(nCh); if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries()); nT = GetListOfTracks(&genParticles,fTrackTypeGen); @@ -690,56 +1017,141 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) vector inputParticlesRec; vector inputParticlesRecRan; - - - // create a random jet within the acceptance - - if(fUseBackgroundCalc){ - Double_t etaMax = 0.9 - fRparam; - Int_t nCone = 0; - Int_t nConeRan = 0; - Double_t pT = 1; // small number - for(int ir = 0;ir < fNRandomCones;ir++){ - Double_t eta = etaMax*2.*(fRandom->Rndm()-0.5); // +- etamax - Double_t phi = fRandom->Rndm()*2.*TMath::Pi(); // 0 - 2pi - // massless jet - 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); - AliAODJet tmpRec (pX,pY,pZ, p); - tmpRec.SetBgEnergy(0,0); // this is use as temporary storage of the summed p_T below - if(rConeArrayRan)new ((*rConeArrayRan)[nConeRan++]) AliAODJet(tmpRec); - if(rConeArray)new ((*rConeArray)[nCone++]) AliAODJet(tmpRec); - } - } - // 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 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(fUseBackgroundCalc&&rConeArray){ - for(int ir = 0;ir < fNRandomCones;ir++){ - AliAODJet *jC = (AliAODJet*)rConeArray->At(ir); - if(jC&&jC->DeltaR(vp)SetBgEnergy(jC->ChargedBgEnergy()+vp->Pt(),0); + 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 sumEff = 0.; + Double_t pT = 0.; + Double_t eff[3] = {0.}; + Int_t cat[3] = {0}; + + Double_t rnd = fRandom->Uniform(1.); + if( fEfficiencyFixed<1. ) { + sumEff = fEfficiencyFixed; + if (fUseDiceEfficiency == 2) { + sumEff = (nCh - fEfficiencyFixed*nGen) / nCh; // rescale eff; fEfficiencyFixed is wrt to nGen, but dicing is fraction of nCh + } + } else { + + pT = vp->Pt(); + 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)); + + //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; + } + } + + 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 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)); @@ -753,67 +1165,19 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) jInpRan.set_user_index(i); inputParticlesRecRan.push_back(jInpRan); vTmpRan.SetPxPyPzE(pX,pY,pZ,p); - - if(fUseBackgroundCalc&&rConeArrayRan){ - for(int ir = 0;ir < fNRandomCones;ir++){ - AliAODJet *jC = (AliAODJet*)rConeArrayRan->At(ir); - if(jC&&jC->DeltaR(&vTmpRan)SetBgEnergy(jC->ChargedBgEnergy()+vTmpRan.Pt(),0); - } - } - } } // 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)); + 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 - if(fUseBackgroundCalc){ - 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 eta = rC->Eta(); - Double_t phi = rC->Phi(); - // massless jet, unit vector - Double_t pT = rC->ChargedBgEnergy(); - if(pT<=0)pT = 0.1; // for almost empty events - 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); - rC->SetPxPyPzE(pX,pY,pZ, p); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - rC = (AliAODJet*)rConeArrayRan->At(ir); - // same wit random - if(rC){ - Double_t eta = rC->Eta(); - Double_t phi = rC->Phi(); - // massless jet, unit vector - Double_t pT = rC->ChargedBgEnergy(); - if(pT<=0)pT = 0.1;// for almost empty events - 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); - rC->SetPxPyPzE(pX,pY,pZ, p); - rC->SetBgEnergy(0,0); - rC->SetEffArea(jetArea,0); - } - } - } - if(inputParticlesRec.size()==0){ if(fDebug)Printf("%s:%d No input particles found, skipping event",(char*)__FILE__,__LINE__); @@ -826,17 +1190,18 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // now create the object that holds info about ghosts - if(!fUseBackgroundCalc&& fNonStdBranch.Length()==0){ + /* + 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 inclusiveJets, sortedJets; fastjet::ClusterSequenceArea clustSeq(inputParticlesRec, jetDef,areaDef); //range where to compute background @@ -848,19 +1213,19 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fastjet::RangeDefinition range(rapMin,rapMax, phiMin, phiMax); + const vector &inclusiveJets = clustSeq.inclusive_jets(); + const vector &sortedJets = sorted_by_pt(inclusiveJets); - inclusiveJets = clustSeq.inclusive_jets(); - sortedJets = sorted_by_pt(inclusiveJets); fh1NJetsRec->Fill(sortedJets.size()); - // loop over all jets an fill information, first on is the leading jet + // loop over all jets an fill information, first on is the leading jet Int_t nRecOver = inclusiveJets.size(); Int_t nRec = inclusiveJets.size(); if(inclusiveJets.size()>0){ AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E()); - Double_t area=clustSeq.area(sortedJets[0]); + Double_t area = clustSeq.area(sortedJets[0]); leadingJet.SetEffArea(area,0); Float_t pt = leadingJet.Pt(); Int_t nAodOutJets = 0; @@ -886,7 +1251,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) Float_t pTback = 0; if(externalBackground){ // carefull has to be filled in a task before - // todo, ReArrange to the botom + // todo, ReArrange to the botom pTback = externalBackground->GetBackground(2)*leadingJet.EffectiveAreaCharged(); } pt = leadingJet.Pt() - pTback; @@ -904,100 +1269,374 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) 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); + } + } + 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&&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); + } + + + 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 constituentsemacs + const vector &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 + + 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){ + 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(partLead); + if(aodT){ + if(aodT->TestFilterBit(fFilterMaskBestPt)){ + aodOutJet->SetTrigger(AliAODJet::kHighTrackPtBest); + } + } + } + } - 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 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); - } + // correlation + Float_t tmpPhi = tmpRec.Phi(); + Float_t tmpEta = tmpRec.Eta(); + if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; + if(j==0){ + fh1PtJetsLeadingRecIn->Fill(tmpPt); + fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); + fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); + fh1NConstLeadingRec->Fill(constituents.size()); + fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); + continue; + } + fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); + fh2JetEtaPt->Fill(tmpEta,tmpPt); + Float_t dPhi = phi - tmpPhi; + if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); + if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); + Float_t dEta = eta - tmpRec.Eta(); + fh2JetsLeadingPhiEta->Fill(dPhi,dEta); + fh2JetsLeadingPhiPt->Fill(dPhi,pt); + if(cenClass>=0){ + fh2JetsLeadingPhiPtC[cenClass]->Fill(dPhi,pt); + 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;iicAt(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 + tmpRecC.SetPtLeading(-1.); + if(fTCARandomConesOut)new ((*fTCARandomConesOut)[nCone++]) AliAODJet(tmpRecC); + if(fTCARandomConesOutRan)new ((*fTCARandomConesOutRan)[nConeRan++]) AliAODJet(tmpRecC); + }// loop over random cones creation - - 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()); - } - - // correlation - Float_t tmpPhi = tmpRec.Phi(); - Float_t tmpEta = tmpRec.Eta(); - if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; - - 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... - - if(j==0){ - fh1PtJetsLeadingRecIn->Fill(tmpPt); - fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta); - fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt); - fh1NConstLeadingRec->Fill(constituents.size()); - fh2NConstLeadingPt->Fill(tmpPt,constituents.size()); - continue; - } - fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta); - fh2JetEtaPt->Fill(tmpEta,tmpPt); - Float_t dPhi = phi - tmpPhi; - if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); - if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); - Float_t dEta = eta - tmpRec.Eta(); - fh2JetsLeadingPhiEta->Fill(dPhi,dEta); - fh2JetsLeadingPhiPt->Fill(dPhi,pt); - fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt); - } - delete recIter; - //background estimates:all bckg jets(0) & wo the 2 hardest(1) - - if(evBkg){ - vector jets2=sortedJets; - if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); - Double_t bkg1=0; - Double_t sigma1=0.; - Double_t meanarea1=0.; - Double_t bkg2=0; - Double_t sigma2=0.; - Double_t meanarea2=0.; - - 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)); + // 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)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(); + 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)fMaxTrackPtInJet)jC->SetTrigger(AliAODJet::kHighTrackPtTriggered); + jC->SetBgEnergy(jC->ChargedBgEnergy()+vTmpRanR.Pt(),0); + if(vTmpRanR.Pt() > jC->GetPtLeading()) jC->SetPtLeading(vTmpRanR.Pt()); + } + } + } + } + }// loop over recparticles + + Float_t jetArea = fRparam*fRparam*TMath::Pi(); + if(fTCARandomConesOut){ + for(int ir = 0;ir < fTCARandomConesOut->GetEntriesFast();ir++){ + // rescale the momentum 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(fAODJetBackgroundOut){ + vector jets2=sortedJets; + if(jets2.size()>2) jets2.erase(jets2.begin(),jets2.begin()+2); + + Double_t bkg1=0; + Double_t sigma1=0.; + Double_t meanarea1=0.; + Double_t bkg2=0; + Double_t sigma2=0.; + Double_t meanarea2=0.; + + Double_t bkg4=0; + Double_t sigma4=0.; + Double_t meanarea4=0.; + + 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(sortedJets, range, true, bkg4, sigma4, meanarea4, true); + fAODJetBackgroundOut->SetBackground(3,bkg4,sigma4,meanarea4); + + ////////////////////// - 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)); + clustSeq.get_median_rho_and_sigma(jets2, range, false, bkg2, sigma2, meanarea2, true); + fAODJetBackgroundOut->SetBackground(1,bkg2,sigma2,meanarea2); + // fh1BiARandomCones[1]->Fill(ptRandomCone-(bkg2*areaRandomCone)); + // fh1BiARandomConesRan[1]->Fill(ptRandomConeRan-(bkg2*areaRandomCone)); + + } + } + + if(fStoreRhoLeadingTrackCorr) { + vector 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); + + //Get leading track in event + AliVParticle *leading = (AliVParticle*)recParticles.At(0); + + fh3CentvsRhoLeadingTrackPt->Fill(cent,bkg2,leading->Pt()); + fh3CentvsSigmaLeadingTrackPt->Fill(cent,sigma2,leading->Pt()); + fh3MultvsRhoLeadingTrackPt->Fill(nCh,bkg2,leading->Pt()); + fh3MultvsSigmaLeadingTrackPt->Fill(nCh,sigma2,leading->Pt()); + + + //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 + 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()); + + fh3CentvsDeltaRhoLeadingTrackPtQ1->Fill(cent,bkg2Q[0]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ2->Fill(cent,bkg2Q[1]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ3->Fill(cent,bkg2Q[2]-bkg2,leading->Pt()); + fh3CentvsDeltaRhoLeadingTrackPtQ4->Fill(cent,bkg2Q[3]-bkg2,leading->Pt()); - } } @@ -1009,77 +1648,79 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // do the same for tracks and jets 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); - while(ptNext()); - if(tmpRec){ - pt = tmpRec->Pt(); - } - } - if(nTrackOver<=0)break; - fh2NRecTracksPt->Fill(ptCut,nTrackOver); - } + 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); + while(ptNext()); + if(tmpRec){ + pt = tmpRec->Pt(); + } + } + if(nTrackOver<=0)break; + fh2NRecTracksPt->Fill(ptCut,nTrackOver); + } - recIter->Reset(); - AliVParticle *leading = (AliVParticle*)recParticles.At(0); - Float_t phi = leading->Phi(); - if(phi<0)phi+=TMath::Pi()*2.; - Float_t eta = leading->Eta(); - pt = leading->Pt(); - while((tmpRec = (AliVParticle*)(recIter->Next()))){ - Float_t tmpPt = tmpRec->Pt(); - Float_t tmpEta = tmpRec->Eta(); - fh1PtTracksRecIn->Fill(tmpPt); - fh2TrackEtaPt->Fill(tmpEta,tmpPt); - if(tmpRec==leading){ - fh1PtTracksLeadingRecIn->Fill(tmpPt); - fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); - continue; - } + recIter->Reset(); + AliVParticle *leading = (AliVParticle*)recParticles.At(0); + Float_t phi = leading->Phi(); + if(phi<0)phi+=TMath::Pi()*2.; + Float_t eta = leading->Eta(); + pt = leading->Pt(); + while((tmpRec = (AliVParticle*)(recIter->Next()))){ + Float_t tmpPt = tmpRec->Pt(); + Float_t tmpEta = tmpRec->Eta(); + fh1PtTracksRecIn->Fill(tmpPt); + fh2TrackEtaPt->Fill(tmpEta,tmpPt); + if(tmpRec==leading){ + fh1PtTracksLeadingRecIn->Fill(tmpPt); + fh2LeadingTrackEtaPt->Fill(tmpEta,tmpPt); + continue; + } // correlation - Float_t tmpPhi = tmpRec->Phi(); + Float_t tmpPhi = tmpRec->Phi(); - if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; - Float_t dPhi = phi - tmpPhi; - if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); - if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); - Float_t dEta = eta - tmpRec->Eta(); - fh2TracksLeadingPhiEta->Fill(dPhi,dEta); - fh2TracksLeadingPhiPt->Fill(dPhi,pt); - fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); - } - delete recIter; - } + if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.; + Float_t dPhi = phi - tmpPhi; + if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi(); + if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi(); + Float_t dEta = eta - tmpRec->Eta(); + fh2TracksLeadingPhiEta->Fill(dPhi,dEta); + fh2TracksLeadingPhiPt->Fill(dPhi,pt); + fh2TracksLeadingPhiPtW->Fill(dPhi,pt,tmpPt); + } + delete recIter; + } - // find the random jets - vector inclusiveJetsRan, sortedJetsRan; - fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); - - inclusiveJetsRan = clustSeqRan.inclusive_jets(); - sortedJetsRan = sorted_by_pt(inclusiveJetsRan); + // find the random jets - // fill the jet information from random track + fastjet::ClusterSequenceArea clustSeqRan(inputParticlesRecRan, jetDef, areaDef); + + // fill the jet information from random track + const vector &inclusiveJetsRan = clustSeqRan.inclusive_jets(); + const vector &sortedJetsRan = sorted_by_pt(inclusiveJetsRan); fh1NJetsRecRan->Fill(sortedJetsRan.size()); - // loop over all jets an fill information, first on is the leading jet + // loop over all jets an fill information, first on is the leading jet + + Int_t nRecOverRan = inclusiveJetsRan.size(); + Int_t nRecRan = inclusiveJetsRan.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(); + 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; - for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ - Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); + Int_t iCount = 0; + TLorentzVector vecarearanb; + + for(int i = 1;i <= fh2NRecJetsPtRan->GetNbinsX();i++){ + Float_t ptCut = fh2NRecJetsPtRan->GetXaxis()->GetBinCenter(i); while(ptFill(tmpPt); // Fill Spectra with constituents - vector constituents = clustSeqRan.constituents(sortedJetsRan[j]); + const vector &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); - Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); - - aodOutJetRan->SetEffArea(arearan,0); } - - - - - for(unsigned int ic = 0; ic < constituents.size();ic++){ - // AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index()); - // fh1PtJetConstRec->Fill(part->Pt()); - if(aodOutJetRan){ - aodOutJetRan->AddTrack(fRef->At(constituents[ic].user_index())); - } - } - - + if(tmpPt>fJetOutputMinPt&&fTCAJetsOutRan){ + aodOutJetRan = new ((*fTCAJetsOutRan)[nAodOutJetsRan++]) AliAODJet(tmpRec); + Double_t arearan=clustSeqRan.area(sortedJetsRan[j]); + 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); + } // correlation Float_t tmpPhi = tmpRec.Phi(); @@ -1156,90 +1788,223 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) } - if(evBkg){ - Double_t bkg3=0.; - Double_t sigma3=0.; - Double_t meanarea3=0.; - clustSeqRan.get_median_rho_and_sigma(sortedJetsRan ,range, false, bkg3, sigma3, meanarea3, true); - evBkg->SetBackground(2,bkg3,sigma3,meanarea3); - // float areaRandomCone = rRandomCone2 *TMath::Pi(); - /* - fh1BiARandomCones[2]->Fill(ptRandomCone-(bkg3*areaRandomCone)); - fh1BiARandomConesRan[2]->Fill(ptRandomConeRan-(bkg3*areaRandomCone)); - */ + 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); + fAODJetBackgroundOut->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 cent = 100; - if(fAOD)cent = fAOD->GetHeader()->GetCentrality(); - Float_t minPt = 9999999; - // hard coded for now ... - // 54.50 44.5 29.5 18.5 for anti-kt rejection 10E-3 - if(cent<10)minPt = 50; - else if(cent<20)minPt = 42; - else if(cent<50)minPt = 28; - else if(cent<80)minPt = 18; - fh1Centrality->Fill(cent); - - float rho = 0; - if(externalBackground)rho = externalBackground->GetBackground(2); - if(jarray&¢<80){ - 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(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); - fh1CentralitySelect->Fill(cent); - aodH->SetFillAOD(kTRUE); - } - } - if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); - PostData(1, fHistList); + // 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(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; + break; + } + } + } + + if(select){ + static AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); + fh1CentralitySelect->Fill(cent); + fh1ZSelect->Fill(zVtx); + aodH->SetFillAOD(kTRUE); + } + } + 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); } void AliAnalysisTaskJetCluster::Terminate(Option_t */*option*/) { -// Terminate analysis -// - if (fDebug > 1) printf("AnalysisJetCluster: Terminate() \n"); + // + // 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){ - if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,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(type==kTrackAOD){ - AliAODEvent *aod = 0; - if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); - else aod = AODEvent(); - if(!aod){ - return iCount; + if(type==kTrackAOD || type==kTrackAODextra || type==kTrackAODextraonly || type==kTrackAODMCextra || type==kTrackAODMCextraonly){ + + if(type!=kTrackAODextraonly && type!=kTrackAODMCextraonly) { + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(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; + } + + // heavy flavor jets + if(fFilterMask==528 && fUseHFcuts){ + Double_t ntpcClus = tr->GetTPCNcls(); + Double_t trPt=tr->Pt(); + TFormula NTPCClsCut("f1NClustersTPCLinearPtDep","70.+30./20.*x"); + + if (trPt <= 20. && (ntpcClus < NTPCClsCut.Eval(trPt))) continue; + else if (trPt > 20. && ntpcClus < 100) continue; + + if (AvoidDoubleCountingHF(aod,tr)) continue; + } + // + + if(fRequireITSRefit){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;} + if (fApplySharedClusterCut) { + Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls()); + if (frac > 0.4) 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()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()Add(tr); - iCount++; + if(type==kTrackAODextra || type==kTrackAODextraonly) { + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + + if(!aod){ + return iCount; + } + TClonesArray *aodExtraTracks = dynamic_cast(aod->FindListObject("aodExtraTracks")); + if(!aodExtraTracks)return iCount; + for(int it =0; itGetEntries(); it++) { + AliVParticle *track = dynamic_cast ((*aodExtraTracks)[it]); + if (!track) continue; + + AliAODTrack *trackAOD = dynamic_cast (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(fRequireITSRefit){if((trackAOD->GetStatus()&AliESDtrack::kITSrefit)==0)continue;} + if (fApplySharedClusterCut) { + Double_t frac = Double_t(trackAOD->GetTPCnclsS()) /Double_t(trackAOD->GetTPCncls()); + if (frac > 0.4) continue; + } + + + if(TMath::Abs(trackAOD->Eta())>fTrackEtaWindow) continue; + if(trackAOD->Pt()Pt()); + list->Add(trackAOD); + iCount++; + } } + + if(type==kTrackAODMCextra || type==kTrackAODMCextraonly) { //embed generator level particles + AliAODEvent *aod = 0; + if(fUseAODTrackInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod){ + return iCount; + } + TClonesArray *aodExtraTracks = dynamic_cast(aod->FindListObject("aodExtraMCparticles")); + if(!aodExtraTracks)return iCount; + for(int it =0; itGetEntries(); it++) { + AliVParticle *track = dynamic_cast ((*aodExtraTracks)[it]); + AliAODMCParticle *partmc = dynamic_cast ((*aodExtraTracks)[it]); + if (!track) { + if(fDebug>10) printf("track %d does not exist\n",it); + continue; + } + + if(!partmc) continue; + if(!partmc->IsPhysicalPrimary())continue; + + if (track->Pt()10) printf("track %d has too low pt %.2f\n",it,track->Pt()); + continue; + } + + /* + AliAODTrack *trackAOD = dynamic_cast((*aodExtraTracks)[it]);//(track); + + if(!trackAOD) { + if(fDebug>10) printf("trackAOD %d does not exist\n",it); + continue; + } + + trackAOD->SetFlags(AliESDtrack::kEmbedded); + trackAOD->SetFilterMap(fFilterMask); + */ + if(fDebug>10) printf("pt extra track %.2f \n", track->Pt()); + + if(TMath::Abs(track->Eta())>fTrackEtaWindow) continue; + if(track->Pt()Add(track); + + iCount++; + } + } + } else if (type == kTrackKineAll||type == kTrackKineCharged){ AliMCEvent* mcEvent = MCEvent(); @@ -1283,28 +2048,255 @@ Int_t AliAnalysisTaskJetCluster::GetListOfTracks(TList *list,Int_t type){ list->Add(part); } else { - if(TMath::Abs(part->Eta())>0.9)continue; + if(TMath::Abs(part->Eta())>fTrackEtaWindow)continue; list->Add(part); } iCount++; } } }// AODMCparticle + else if (type == kTrackAODMCHF){ + + AliAODEvent *aod = 0; + if(fUseAODMCInput)aod = dynamic_cast(InputEvent()); + else aod = AODEvent(); + if(!aod)return iCount; + TClonesArray *tca = dynamic_cast(aod->FindListObject(AliAODMCParticle::StdBranchName())); + if(!tca)return iCount; + for(int it = 0;it < tca->GetEntriesFast();++it){ + AliAODMCParticle *part = (AliAODMCParticle*)(tca->At(it)); + if(!part) continue; + int partpdg= part->PdgCode(); + if(!part->IsPhysicalPrimary() && !IsBMeson(partpdg) && !IsDMeson(partpdg) )continue; + + if (IsBMeson(partpdg) || IsDMeson(partpdg)) { + iCount+= AddDaughters( list , part,tca); + } + else { + + if(part->Pt()Eta())>fTrackEtaWindow) continue; + if(part->Charge()==0) continue; + + if((part->Pt()>=fTrackPtCut) && (TMath::Abs(part->Eta())<=fTrackEtaWindow) && (part->Charge()!=0))list->Add(part); + iCount++; + } + } + } + list->Sort(); return iCount; } +Int_t AliAnalysisTaskJetCluster::AddDaughters(TList * list, AliAODMCParticle *part, TClonesArray * tca){ + Int_t count=0; + Int_t nDaugthers = part->GetNDaughters(); + for(Int_t i=0;i< nDaugthers;++i){ + AliAODMCParticle *partdaughter = (AliAODMCParticle*)(tca->At(i)); + if(!partdaughter) continue; + if(partdaughter->Pt()Eta())>fTrackEtaWindow)continue; + if(partdaughter->Charge()==0)continue; + + if(!IsDMeson(partdaughter->PdgCode()) && !IsBMeson(partdaughter->PdgCode()) ){ + list->Add(partdaughter); + count++; + } + else count+=AddDaughters(list,part,tca); + } +return count; +} + + +Bool_t AliAnalysisTaskJetCluster::AvoidDoubleCountingHF(AliAODEvent *aod, AliAODTrack *tr1){ + + if(!(tr1->TestFilterBit(BIT(9)))) return kFALSE; + + Int_t idtr1 = tr1->GetID(); + + for(int jt = 0;jt < aod->GetNumberOfTracks();++jt){ + + const AliAODTrack *tr2 = aod->GetTrack(jt); + Int_t idtr2 = tr2->GetID(); + + if (!(tr2->TestFilterBit(BIT(4)))) continue; + if (idtr1==(idtr2+1)*-1.) return kTRUE; + + } + return kFALSE; +} +void AliAnalysisTaskJetCluster::LoadTrPtResolutionRootFileFromOADB() { + + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + + 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() { + + if (!gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + + 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 &inputParticles){ + Int_t AliAnalysisTaskJetCluster::AddParticlesFastJet(TList &particles,vector &inputParticles){ for(int i = 0; i < particles.GetEntries(); i++){ - AliVParticle *vp = (AliVParticle*)particles.At(i); - // Carefull energy is not well determined in real data, should not matter for p_T scheme? - fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); - jInp.set_user_index(i); - inputParticles.push_back(jInp); + AliVParticle *vp = (AliVParticle*)particles.At(i); + // Carefull energy is not well determined in real data, should not matter for p_T scheme? + fastjet::PseudoJet jInp(vp->Px(),vp->Py(),vp->Pz(),vp->E()); + jInp.set_user_index(i); + inputParticles.push_back(jInp); } return 0; -} + } */ + + +bool AliAnalysisTaskJetCluster::IsBMeson(int pc){ + int bPdG[] = {511,521,10511,10521,513,523,10513,10523,20513,20523,20513,20523,515,525,531, + 10531,533,10533,20533,535,541,10541,543,10543,20543,545,551,10551,100551, + 110551,200551,210551,553,10553,20553,30553,100553,110553,120553,130553,200553,210553,220553, + 300553,9000533,9010553,555,10555,20555,100555,110555,120555,200555,557,100557}; + for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true; +return false; +} +bool AliAnalysisTaskJetCluster::IsDMeson(int pc){ + int bPdG[] = {411,421,10411,10421,413,423,10413,10423,20431,20423,415, + 425,431,10431,433,10433,20433,435,441,10441,100441,443,10443,20443, + 100443,30443,9000443,9010443,9020443,445,100445}; + for(int i=0;i< (int)(sizeof(bPdG)/sizeof(int));++i) if(abs(pc) == bPdG[i]) return true; +return false; +} +