From: kleinb Date: Thu, 22 Dec 2011 12:30:15 +0000 (+0000) Subject: fragmentation function correction task added (O. Busch) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=39e2b0572f53eb082d43a612ca02a21ba41b3eeb;p=u%2Fmrichter%2FAliRoot.git fragmentation function correction task added (O. Busch) --- diff --git a/PWG4/CMakelibPWG4JetTasks.pkg b/PWG4/CMakelibPWG4JetTasks.pkg index 010c5877bfb..aaafb5feede 100644 --- a/PWG4/CMakelibPWG4JetTasks.pkg +++ b/PWG4/CMakelibPWG4JetTasks.pkg @@ -43,7 +43,8 @@ set ( SRCS JetTasks/AliPWG4HighPtSpectra.cxx JetTasks/AliPWG4CosmicCandidates.cxx JetTasks/AliAnalysisTaskJetChem.cxx - JetTasks/AliAnalysisTaskFragmentationFunction.cxx + JetTasks/AliAnalysisTaskFragmentationFunction.cxx + JetTasks/AliFragmentationFunctionCorrections.cxx JetTasks/AliAnalysisTaskMinijet.cxx JetTasks/AliAnalysisTaskQGSep.cxx JetTasks/AliAnalysisTaskJetCore.cxx diff --git a/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx b/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx index cfe745784b7..493ed7e739c 100644 --- a/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx @@ -65,6 +65,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fBranchRecBackJets("") ,fBranchRecBckgClusters("") ,fBranchGenJets("") + ,fBranchEmbeddedJets("") ,fTrackTypeGen(0) ,fJetTypeGen(0) ,fJetTypeRecEff(0) @@ -78,6 +79,11 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fTrackEtaMax(0) ,fTrackPhiMin(0) ,fTrackPhiMax(0) + ,fUseExtraTracks(0) + ,fUseExtraTracksBgr(0) + ,fCutFractionPtEmbedded(0) + ,fUseEmbeddedJetAxis(0) + ,fUseEmbeddedJetPt(0) ,fJetPtCut(0) ,fJetEtaMin(0) ,fJetEtaMax(0) @@ -110,6 +116,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fJetsRecCuts(0) ,fJetsGen(0) ,fJetsRecEff(0) + ,fJetsEmbedded(0) ,fBckgJetsRec(0) ,fBckgJetsRecCuts(0) ,fBckgJetsGen(0) @@ -251,6 +258,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fh1nRecJetsCuts(0) ,fh1nGenJets(0) ,fh1nRecEffJets(0) + ,fh1nEmbeddedJets(0) ,fh1nRecBckgJetsCuts(0) ,fh1nGenBckgJets(0) ,fh2PtRecVsGenPrim(0) @@ -261,7 +269,8 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fFFHistosRecEffGen(0) ,fFFHistosRecEffRec(0) ,fFFHistosSecRec(0) - ,fhnResponseSinglePt(0) + ,fhnResponseSinglePt(0) + ,fh2SingleInvPtRecMnGenVsPtGen(0) ,fhnResponseJetTrackPt(0) ,fhnResponseJetZ(0) ,fhnResponseJetXi(0) @@ -276,6 +285,15 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction() ,fh1Out3JetsMult(0) ,fh1MedianClustersMult(0) ,fh1OutClustersMult(0) + ,fh1FractionPtEmbedded(0) + ,fh1IndexEmbedded(0) + ,fh2DeltaPtVsJetPtEmbedded(0) + ,fh2DeltaPtVsRecJetPtEmbedded(0) + ,fh1DeltaREmbedded(0) + ,fh2ptVsDistNN_pt50_rec(0) + ,fh2ptVsDistNN_pt50_nonRec(0) + ,fh2ptVsDistNN_pt10_rec(0) + ,fh2ptVsDistNN_pt10_nonRec(0) ,fQABckgHisto0RecCuts(0) ,fQABckgHisto0Gen(0) ,fQABckgHisto1RecCuts(0) @@ -349,6 +367,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fBranchRecBackJets("") ,fBranchRecBckgClusters("") ,fBranchGenJets("") + ,fBranchEmbeddedJets("") ,fTrackTypeGen(0) ,fJetTypeGen(0) ,fJetTypeRecEff(0) @@ -362,6 +381,11 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fTrackEtaMax(0) ,fTrackPhiMin(0) ,fTrackPhiMax(0) + ,fUseExtraTracks(0) + ,fUseExtraTracksBgr(0) + ,fCutFractionPtEmbedded(0) + ,fUseEmbeddedJetAxis(0) + ,fUseEmbeddedJetPt(0) ,fJetPtCut(0) ,fJetEtaMin(0) ,fJetEtaMax(0) @@ -394,6 +418,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fJetsRecCuts(0) ,fJetsGen(0) ,fJetsRecEff(0) + ,fJetsEmbedded(0) ,fBckgJetsRec(0) ,fBckgJetsRecCuts(0) ,fBckgJetsGen(0) @@ -535,6 +560,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fh1nRecJetsCuts(0) ,fh1nGenJets(0) ,fh1nRecEffJets(0) + ,fh1nEmbeddedJets(0) ,fh1nRecBckgJetsCuts(0) ,fh1nGenBckgJets(0) ,fh2PtRecVsGenPrim(0) @@ -546,6 +572,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fFFHistosRecEffRec(0) ,fFFHistosSecRec(0) ,fhnResponseSinglePt(0) + ,fh2SingleInvPtRecMnGenVsPtGen(0) ,fhnResponseJetTrackPt(0) ,fhnResponseJetZ(0) ,fhnResponseJetXi(0) @@ -560,6 +587,15 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fh1Out3JetsMult(0) ,fh1MedianClustersMult(0) ,fh1OutClustersMult(0) + ,fh1FractionPtEmbedded(0) + ,fh1IndexEmbedded(0) + ,fh2DeltaPtVsJetPtEmbedded(0) + ,fh2DeltaPtVsRecJetPtEmbedded(0) + ,fh1DeltaREmbedded(0) + ,fh2ptVsDistNN_pt50_rec(0) + ,fh2ptVsDistNN_pt50_nonRec(0) + ,fh2ptVsDistNN_pt10_rec(0) + ,fh2ptVsDistNN_pt10_nonRec(0) ,fQABckgHisto0RecCuts(0) ,fQABckgHisto0Gen(0) ,fQABckgHisto1RecCuts(0) @@ -637,6 +673,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fBranchRecBackJets(copy.fBranchRecBackJets) ,fBranchRecBckgClusters(copy.fBranchRecBckgClusters) ,fBranchGenJets(copy.fBranchGenJets) + ,fBranchEmbeddedJets(copy.fBranchEmbeddedJets) ,fTrackTypeGen(copy.fTrackTypeGen) ,fJetTypeGen(copy.fJetTypeGen) ,fJetTypeRecEff(copy.fJetTypeRecEff) @@ -650,6 +687,11 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fTrackEtaMax(copy.fTrackEtaMax) ,fTrackPhiMin(copy.fTrackPhiMin) ,fTrackPhiMax(copy.fTrackPhiMax) + ,fUseExtraTracks(copy.fUseExtraTracks) + ,fUseExtraTracksBgr(copy.fUseExtraTracksBgr) + ,fCutFractionPtEmbedded(copy.fCutFractionPtEmbedded) + ,fUseEmbeddedJetAxis(copy.fUseEmbeddedJetAxis) + ,fUseEmbeddedJetPt(copy.fUseEmbeddedJetPt) ,fJetPtCut(copy.fJetPtCut) ,fJetEtaMin(copy.fJetEtaMin) ,fJetEtaMax(copy.fJetEtaMax) @@ -682,6 +724,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fJetsRecCuts(copy.fJetsRecCuts) ,fJetsGen(copy.fJetsGen) ,fJetsRecEff(copy.fJetsRecEff) + ,fJetsEmbedded(copy.fJetsEmbedded) ,fBckgJetsRec(copy.fBckgJetsRec) ,fBckgJetsRecCuts(copy.fBckgJetsRecCuts) ,fBckgJetsGen(copy.fBckgJetsGen) @@ -823,6 +866,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fh1nRecJetsCuts(copy.fh1nRecJetsCuts) ,fh1nGenJets(copy.fh1nGenJets) ,fh1nRecEffJets(copy.fh1nRecEffJets) + ,fh1nEmbeddedJets(copy.fh1nEmbeddedJets) ,fh1nRecBckgJetsCuts(copy.fh1nRecBckgJetsCuts) ,fh1nGenBckgJets(copy.fh1nGenBckgJets) ,fh2PtRecVsGenPrim(copy.fh2PtRecVsGenPrim) @@ -834,6 +878,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fFFHistosRecEffRec(copy.fFFHistosRecEffRec) ,fFFHistosSecRec(copy.fFFHistosSecRec) ,fhnResponseSinglePt(copy.fhnResponseSinglePt) + ,fh2SingleInvPtRecMnGenVsPtGen(fh2SingleInvPtRecMnGenVsPtGen) ,fhnResponseJetTrackPt(copy.fhnResponseJetTrackPt) ,fhnResponseJetZ(copy.fhnResponseJetZ) ,fhnResponseJetXi(copy.fhnResponseJetXi) @@ -848,6 +893,15 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const ,fh1Out3JetsMult(copy.fh1Out3JetsMult) ,fh1MedianClustersMult(copy.fh1MedianClustersMult) ,fh1OutClustersMult(copy.fh1OutClustersMult) + ,fh1FractionPtEmbedded(copy.fh1FractionPtEmbedded) + ,fh1IndexEmbedded(copy.fh1IndexEmbedded) + ,fh2DeltaPtVsJetPtEmbedded(copy.fh2DeltaPtVsJetPtEmbedded) + ,fh2DeltaPtVsRecJetPtEmbedded(copy.fh2DeltaPtVsRecJetPtEmbedded) + ,fh1DeltaREmbedded(copy.fh1DeltaREmbedded) + ,fh2ptVsDistNN_pt50_rec(copy.fh2ptVsDistNN_pt50_rec) + ,fh2ptVsDistNN_pt50_nonRec(copy.fh2ptVsDistNN_pt50_nonRec) + ,fh2ptVsDistNN_pt10_rec(copy.fh2ptVsDistNN_pt10_rec) + ,fh2ptVsDistNN_pt10_nonRec(copy.fh2ptVsDistNN_pt10_nonRec) ,fQABckgHisto0RecCuts(copy.fQABckgHisto0RecCuts) ,fQABckgHisto0Gen(copy.fQABckgHisto0Gen) ,fQABckgHisto1RecCuts(copy.fQABckgHisto1RecCuts) @@ -926,6 +980,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fBranchRecBackJets = o.fBranchRecBackJets; fBranchRecBckgClusters = o.fBranchRecBckgClusters; fBranchGenJets = o.fBranchGenJets; + fBranchEmbeddedJets = o.fBranchEmbeddedJets; fTrackTypeGen = o.fTrackTypeGen; fJetTypeGen = o.fJetTypeGen; fJetTypeRecEff = o.fJetTypeRecEff; @@ -939,6 +994,11 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fTrackEtaMax = o.fTrackEtaMax; fTrackPhiMin = o.fTrackPhiMin; fTrackPhiMax = o.fTrackPhiMax; + fUseExtraTracks = o.fUseExtraTracks; + fUseExtraTracksBgr = o.fUseExtraTracksBgr; + fCutFractionPtEmbedded = o.fCutFractionPtEmbedded; + fUseEmbeddedJetAxis = o.fUseEmbeddedJetAxis; + fUseEmbeddedJetPt = o.fUseEmbeddedJetPt; fJetPtCut = o.fJetPtCut; fJetEtaMin = o.fJetEtaMin; fJetEtaMax = o.fJetEtaMax; @@ -976,6 +1036,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fJetsRecCuts = o.fJetsRecCuts; fJetsGen = o.fJetsGen; fJetsRecEff = o.fJetsRecEff; + fJetsEmbedded = o.fJetsEmbedded; fBckgJetsRec = o.fBckgJetsRec; fBckgJetsRecCuts = o.fBckgJetsRecCuts; fBckgJetsGen = o.fBckgJetsGen; @@ -1117,6 +1178,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fh1nRecJetsCuts = o.fh1nRecJetsCuts; fh1nGenJets = o.fh1nGenJets; fh1nRecEffJets = o.fh1nRecEffJets; + fh1nEmbeddedJets = o.fh1nEmbeddedJets; fh2PtRecVsGenPrim = o.fh2PtRecVsGenPrim; fh2PtRecVsGenSec = o.fh2PtRecVsGenSec; fQATrackHistosRecEffGen = o.fQATrackHistosRecEffGen; @@ -1126,6 +1188,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fFFHistosRecEffRec = o.fFFHistosRecEffRec; fFFHistosSecRec = o.fFFHistosSecRec; fhnResponseSinglePt = o.fhnResponseSinglePt; + fh2SingleInvPtRecMnGenVsPtGen = o.fh2SingleInvPtRecMnGenVsPtGen; fhnResponseJetTrackPt = o.fhnResponseJetTrackPt; fhnResponseJetZ = o.fhnResponseJetZ; fhnResponseJetXi = o.fhnResponseJetXi; @@ -1140,6 +1203,15 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper fh1Out3JetsMult = o.fh1Out3JetsMult; fh1MedianClustersMult = o.fh1MedianClustersMult; fh1OutClustersMult = o.fh1OutClustersMult; + fh1FractionPtEmbedded = o.fh1FractionPtEmbedded; + fh1IndexEmbedded = o.fh1IndexEmbedded; + fh2DeltaPtVsJetPtEmbedded = o.fh2DeltaPtVsJetPtEmbedded; + fh2DeltaPtVsRecJetPtEmbedded = o.fh2DeltaPtVsRecJetPtEmbedded; + fh1DeltaREmbedded = o.fh1DeltaREmbedded; + fh2ptVsDistNN_pt50_rec = o.fh2ptVsDistNN_pt50_rec; + fh2ptVsDistNN_pt50_nonRec = o.fh2ptVsDistNN_pt50_nonRec; + fh2ptVsDistNN_pt10_rec = o.fh2ptVsDistNN_pt10_rec; + fh2ptVsDistNN_pt10_nonRec = o.fh2ptVsDistNN_pt10_nonRec; fQABckgHisto0RecCuts = o.fQABckgHisto0RecCuts; fQABckgHisto0Gen = o.fQABckgHisto0Gen; fQABckgHisto1RecCuts = o.fQABckgHisto1RecCuts; @@ -1208,6 +1280,8 @@ AliAnalysisTaskFragmentationFunction::~AliAnalysisTaskFragmentationFunction() if(fJetsRecCuts) delete fJetsRecCuts; if(fJetsGen) delete fJetsGen; if(fJetsRecEff) delete fJetsRecEff; + if(fJetsEmbedded) delete fJetsEmbedded; + if(fBckgMode && (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters || fBckgType[0]==kBckgClustersOutLeading || fBckgType[1]==kBckgClustersOutLeading || fBckgType[2]==kBckgClustersOutLeading || @@ -2352,6 +2426,9 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() fJetsRecEff = new TList(); fJetsRecEff->SetOwner(kFALSE); + fJetsEmbedded = new TList(); + fJetsEmbedded->SetOwner(kFALSE); + // fJetsKine = new TList(); // fJetsKine->SetOwner(kTRUE); // delete AOD jets using mom from Kine Tree via TList::Clear() @@ -2407,9 +2484,12 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5); fh1nGenJets = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5); fh1nRecEffJets = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5); + fh1nEmbeddedJets = new TH1F("fh1nEmbeddedJets","embedded jets per event",10,-0.5,9.5); + fh2PtRecVsGenPrim = new TH2F("fh2PtRecVsGenPrim","rec vs gen pt",fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax,fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax); fh2PtRecVsGenSec = new TH2F("fh2PtRecVsGenSec","rec vs gen pt",fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax,fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax); + // Background if(fBckgMode) { if(fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters || @@ -2443,6 +2523,23 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() fh1OutClustersMult = new TH1F("fh1OutClustersMult","Background multiplicity - clusters outside leading jet",3000,0.,3000.); } + // embedding + if(fBranchEmbeddedJets.Length()){ + fh1FractionPtEmbedded = new TH1F("fh1FractionPtEmbedded","",200,0,2); + fh1IndexEmbedded = new TH1F("fh1IndexEmbedded","",11,-1,10); + fh2DeltaPtVsJetPtEmbedded = new TH2F("fh2DeltaPtVsJetPtEmbedded","",250,0,250,200,-100,100); + fh2DeltaPtVsRecJetPtEmbedded = new TH2F("fh2DeltaPtVsRecJetPtEmbedded","",250,0,250,200,-100,100); + fh1DeltaREmbedded = new TH1F("fh1DeltaREmbedded","",50,0,0.5); + fh1nEmbeddedJets = new TH1F("fh1nEmbeddedJets","embedded jets per event",10,-0.5,9.5); + } + + if(fEffMode){ + fh2ptVsDistNN_pt50_rec = new TH2F("fh2ptVsDistNN_pt50_rec","",200,0,0.2,500,0.,100); + fh2ptVsDistNN_pt50_nonRec = new TH2F("fh2ptVsDistNN_pt50_nonRec","",200,0,0.2,500,0.,100); + fh2ptVsDistNN_pt10_rec = new TH2F("fh2ptVsDistNN_pt10_rec","",200,0,0.2,500,0.,100); + fh2ptVsDistNN_pt10_nonRec = new TH2F("fh2ptVsDistNN_pt10_nonRec","",200,0,0.2,500,0.,100); + } + if(fQAMode){ if(fQAMode&1){ // track QA fQATrackHistosRec = new AliFragFuncQATrackHistos("Rec", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, @@ -2635,6 +2732,9 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() nBinsResponseSinglePt,binMinResponseSinglePt,binMaxResponseSinglePt); AliAnalysisTaskFragmentationFunction::SetProperties(fhnResponseSinglePt,2,labelsResponseSinglePt); + + // TH2F inv pt diff + fh2SingleInvPtRecMnGenVsPtGen = new TH2F("fh2SingleInvPtRecMnGenVsPtGen","",fQATrackNBinsPt,fQATrackPtMin,fQATrackPtMax,200,-1,1); } if(fFFMode){ fFFHistosRecEffGen = new AliFragFuncHistos("RecEffGen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, @@ -3167,6 +3267,23 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() fCommonHistList->Add(fh1OutClustersMult); } + + if(fBranchEmbeddedJets.Length()){ + fCommonHistList->Add(fh1FractionPtEmbedded); + fCommonHistList->Add(fh1IndexEmbedded); + fCommonHistList->Add(fh2DeltaPtVsJetPtEmbedded); + fCommonHistList->Add(fh2DeltaPtVsRecJetPtEmbedded); + fCommonHistList->Add(fh1DeltaREmbedded); + fCommonHistList->Add(fh1nEmbeddedJets); + } + + if(fEffMode){ + fCommonHistList->Add(fh2ptVsDistNN_pt50_rec); + fCommonHistList->Add(fh2ptVsDistNN_pt50_nonRec); + fCommonHistList->Add(fh2ptVsDistNN_pt10_rec); + fCommonHistList->Add(fh2ptVsDistNN_pt10_nonRec); + } + // QA if(fQAMode){ if(fQAMode&1){ // track QA @@ -3260,6 +3377,7 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() fQATrackHistosRecEffRec->AddToOutput(fCommonHistList); fQATrackHistosSecRec->AddToOutput(fCommonHistList); fCommonHistList->Add(fhnResponseSinglePt); + fCommonHistList->Add(fh2SingleInvPtRecMnGenVsPtGen); } if(fFFMode){ fFFHistosRecEffGen->AddToOutput(fCommonHistList); @@ -3286,6 +3404,8 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects() } TH1::AddDirectory(oldStatus); + + PostData(1, fCommonHistList); } //_______________________________________________ @@ -3387,6 +3507,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) } else { cl = AliAnalysisHelperJetTasks::EventClass(); + if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // retrieve value 'by hand' } if(cl!=fEventClass){ @@ -3403,6 +3524,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) Int_t nTracksPrim = primVtx->GetNContributors(); fh1VertexNContributors->Fill(nTracksPrim); + if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim); if(!nTracksPrim){ if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); @@ -3440,34 +3562,37 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) Double_t nTrials = 1; // trials for MC trigger weight for real data if(fMCEvent){ - AliGenEventHeader* genHeader = fMCEvent->GenEventHeader(); - if(genHeader){ - AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(genHeader); - AliGenHijingEventHeader* hijingGenHeader = 0x0; - - if(pythiaGenHeader){ - if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__); - nTrials = pythiaGenHeader->Trials(); - ptHard = pythiaGenHeader->GetPtHard(); - - fh1PtHard->Fill(ptHard); - fh1PtHardTrials->Fill(ptHard,nTrials); - - - } else { // no pythia, hijing? - - if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__); - - hijingGenHeader = dynamic_cast(genHeader); - if(!hijingGenHeader){ - Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__); - } else { - if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__); - } - } + AliGenEventHeader* genHeader = fMCEvent->GenEventHeader(); - fh1Trials->Fill("#sum{ntrials}",fAvgTrials); + if(genHeader){ + + AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast(genHeader); + AliGenHijingEventHeader* hijingGenHeader = 0x0; + + if(pythiaGenHeader){ + if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__); + nTrials = pythiaGenHeader->Trials(); + ptHard = pythiaGenHeader->GetPtHard(); + + fh1PtHard->Fill(ptHard); + fh1PtHardTrials->Fill(ptHard,nTrials); + + + } else { // no pythia, hijing? + + if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__); + + hijingGenHeader = dynamic_cast(genHeader); + if(!hijingGenHeader){ + Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__); + } else { + if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__); + } + } + + fh1Trials->Fill("#sum{ntrials}",fAvgTrials); } + } //___ fetch jets __________________________________________________________________________ @@ -3486,10 +3611,12 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) fh1nRecJetsCuts->Fill(nRecJetsCuts); if(fJetTypeGen==kJetsKine || fJetTypeGen == kJetsKineAcceptance) fJetsGen->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear() + Int_t nJGen = GetListOfJets(fJetsGen, fJetTypeGen); Int_t nGenJets = 0; if(nJGen>=0) nGenJets = fJetsGen->GetEntries(); if(fDebug>2)Printf("%s:%d Selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets); + if(nJGen != nGenJets) Printf("%s:%d Mismatch selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets); fh1nGenJets->Fill(nGenJets); @@ -3502,6 +3629,34 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) if(nJRecEff != nRecEffJets) Printf("%s:%d Mismatch selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets); fh1nRecEffJets->Fill(nRecEffJets); + + Int_t nEmbeddedJets = 0; + TArrayI iEmbeddedMatchIndex; + TArrayF fEmbeddedPtFraction; + + + if(fBranchEmbeddedJets.Length()){ + Int_t nJEmbedded = GetListOfJets(fJetsEmbedded, kJetsEmbedded); + if(nJEmbedded>=0) nEmbeddedJets = fJetsEmbedded->GetEntries(); + if(fDebug>2)Printf("%s:%d Selected Embedded jets: %d %d",(char*)__FILE__,__LINE__,nJEmbedded,nEmbeddedJets); + if(nJEmbedded != nEmbeddedJets) Printf("%s:%d Mismatch Selected Embedded Jets: %d %d",(char*)__FILE__,__LINE__,nJEmbedded,nEmbeddedJets); + fh1nEmbeddedJets->Fill(nEmbeddedJets); + + Float_t maxDist = 0.3; + + iEmbeddedMatchIndex.Set(nEmbeddedJets); + fEmbeddedPtFraction.Set(nEmbeddedJets); + + iEmbeddedMatchIndex.Reset(-1); + fEmbeddedPtFraction.Reset(0); + + AliAnalysisHelperJetTasks::GetJetMatching(fJetsEmbedded, nEmbeddedJets, + fJetsRecCuts, nRecJetsCuts, + iEmbeddedMatchIndex, fEmbeddedPtFraction, + fDebug, maxDist); + + } + //____ fetch background jets ___________________________________________________ if(fBckgMode && (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters || @@ -3535,14 +3690,22 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) //____ fetch particles __________________________________________________________ - Int_t nT = GetListOfTracks(fTracksRec, kTrackAOD); + Int_t nT; + if(fUseExtraTracks == 1) nT = GetListOfTracks(fTracksRec, kTrackAODExtra); + else if(fUseExtraTracks == -1) nT = GetListOfTracks(fTracksRec, kTrackAODExtraonly); + else nT = GetListOfTracks(fTracksRec, kTrackAOD); + Int_t nRecPart = 0; if(nT>=0) nRecPart = fTracksRec->GetEntries(); if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart); if(nRecPart != nT) Printf("%s:%d Mismatch selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart); - Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts); + Int_t nTCuts; + if(fUseExtraTracks == 1) nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODExtraCuts); + else if(fUseExtraTracks == -1) nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODExtraonlyCuts); + else nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts); + Int_t nRecPartCuts = 0; if(nTCuts>=0) nRecPartCuts = fTracksRecCuts->GetEntries(); if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts); @@ -3567,10 +3730,6 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) AliVParticle *part = dynamic_cast(fTracksRec->At(it)); if(part)fQATrackHistosRec->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt()); } - for(Int_t it=0; it(fTracksRecCuts->At(it)); - if(part)fQATrackHistosRecCuts->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt()); - } for(Int_t it=0; it(fTracksGen->At(it)); if(part)fQATrackHistosGen->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt()); @@ -3579,17 +3738,17 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) // fill DCA to prim vertex for(Int_t it=0; it(fTracksRecCuts->At(it)); - + if(!aodtr) continue; if(!primVtx) continue; Double_t Bfield = fAOD->GetMagneticField(); Double_t dz[2]; Double_t cov[3]; - - AliAODTrack tmp(*aodtr); + + AliAODTrack tmp(*aodtr); tmp.PropagateToDCA(primVtx, Bfield, 5., dz, cov); - + Double_t dcaXY = dz[0]; Double_t dcaZ = dz[1]; @@ -3598,9 +3757,6 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) } } - - - // loop over jets if(fQAMode&2){ @@ -3621,102 +3777,147 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) if(fQAMode&2) fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() ); + + Double_t ptFractionEmbedded = 0; + AliAODJet* embeddedJet = 0; + + if(fBranchEmbeddedJets.Length()){ // find embedded jet + + Int_t indexEmbedded = -1; + for(Int_t i=0; iFill(indexEmbedded); + fh1FractionPtEmbedded->Fill(ptFractionEmbedded); + + if(indexEmbedded>-1){ + + embeddedJet = dynamic_cast(fJetsEmbedded->At(indexEmbedded)); + if(!embeddedJet) continue; + + Double_t deltaPt = jet->Pt() - embeddedJet->Pt(); + Double_t deltaR = jet->DeltaR((AliVParticle*) (embeddedJet)); + + fh2DeltaPtVsJetPtEmbedded->Fill(embeddedJet->Pt(),deltaPt); + fh2DeltaPtVsRecJetPtEmbedded->Fill(jet->Pt(),deltaPt); + fh1DeltaREmbedded->Fill(deltaR); + } + } + + // get tracks in jet TList* jettracklist = new TList(); Double_t sumPt = 0.; Float_t leadTrackPt = 0.; TLorentzVector* leadTrackV = new TLorentzVector(); - + if(GetFFRadius()<=0){ GetJetTracksTrackrefs(jettracklist, jet); } else { - GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt); - } - - for(Int_t it=0; itGetSize(); ++it){ - - AliVParticle* trackVP = dynamic_cast(jettracklist->At(it)); - if(!trackVP)continue; - TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); - - Float_t jetPt = jet->Pt(); - Float_t trackPt = trackV->Pt(); - - Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE; - - if(fFFMode) fFFHistosRecCuts->FillFF( trackPt, jetPt, incrementJetPt); - if(fIJMode) fIJHistosRecCuts->FillIntraJet( trackV, jet->MomentumVector() ); - - if(it==0){ // leading track - leadTrackPt = trackPt; - leadTrackV->SetPxPyPzE(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); - - if(fFFMode) fFFHistosRecLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE); - if(fIJMode) fIJHistosRecLeadingTrack->FillIntraJet( leadTrackV, jet->MomentumVector() ); + if(fUseEmbeddedJetAxis){ + if(embeddedJet) GetJetTracksPointing(fTracksRecCuts, jettracklist, embeddedJet, GetFFRadius(), sumPt); } - if(fFFMode) fFFHistosRecLeading->FillFF( trackPt, leadTrackPt , incrementJetPt); - if(fIJMode) fIJHistosRecLeading->FillIntraJet( trackV, leadTrackV ); - - delete trackV; + else GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt); } - // ff and ij for background study - if(fBckgMode){ - if(fBckgType[0]!=-1) - FillBckgHistos(fBckgType[0], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, - fFFBckgHisto0RecCuts, fFFBckgHisto0RecLeading, - fIJBckgHisto0RecCuts, fIJBckgHisto0RecLeading, - fQABckgHisto0RecCuts); - if(fBckgType[1]!=-1) - FillBckgHistos(fBckgType[1], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, - fFFBckgHisto1RecCuts, fFFBckgHisto1RecLeading, - fIJBckgHisto1RecCuts, fIJBckgHisto1RecLeading, - fQABckgHisto1RecCuts); - if(fBckgType[2]!=-1) - FillBckgHistos(fBckgType[2], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, - fFFBckgHisto2RecCuts, fFFBckgHisto2RecLeading, - fIJBckgHisto2RecCuts, fIJBckgHisto2RecLeading, - fQABckgHisto2RecCuts); - if(fBckgType[3]!=-1) - FillBckgHistos(fBckgType[3], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, - fFFBckgHisto3RecCuts, fFFBckgHisto3RecLeading, - fIJBckgHisto3RecCuts, fIJBckgHisto3RecLeading, - fQABckgHisto3RecCuts); - if(fBckgType[4]!=-1) - FillBckgHistos(fBckgType[4], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, - fFFBckgHisto4RecCuts, fFFBckgHisto4RecLeading, - fIJBckgHisto4RecCuts, fIJBckgHisto4RecLeading, - fQABckgHisto4RecCuts); - } // end if(fBckgMode) - - delete leadTrackV; - delete jettracklist; - // phi correlation - if(fPhiCorrMode){ - for(Int_t it=0; itAt(it)); + if(ptFractionEmbedded>=fCutFractionPtEmbedded){ // if no embedding: ptFraction = cutFraction = 0 + + for(Int_t it=0; itGetSize(); ++it){ + + AliVParticle* trackVP = dynamic_cast(jettracklist->At(it)); + if(!trackVP)continue; + TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); - Float_t partEta = part->Eta(); - Float_t partPhi = part->Phi(); - Float_t partPt = part->Pt(); + Float_t jetPt = jet->Pt(); + if(fUseEmbeddedJetPt){ + if(embeddedJet) jetPt = embeddedJet->Pt(); + else jetPt = 0; + } + Float_t trackPt = trackV->Pt(); + + + Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE; - fPhiCorrHistosJetArea->FillTrackQA( partEta, - TVector2::Phi_mpi_pi( jet->Phi() - partPhi ), - partPt, - kTRUE); + if(fFFMode) fFFHistosRecCuts->FillFF( trackPt, jetPt, incrementJetPt); + if(fIJMode) fIJHistosRecCuts->FillIntraJet( trackV, jet->MomentumVector() ); - fPhiCorrHistosTransverseArea->FillTrackQA( partEta, - TVector2::Phi_mpi_pi( jet->Phi() - partPhi + TMath::Pi()/2), - partPt, - kTRUE); + if(it==0){ // leading track + leadTrackPt = trackPt; + leadTrackV->SetPxPyPzE(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); + + if(fFFMode) fFFHistosRecLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE); + if(fIJMode) fIJHistosRecLeadingTrack->FillIntraJet( leadTrackV, jet->MomentumVector() ); + } + if(fFFMode) fFFHistosRecLeading->FillFF( trackPt, leadTrackPt , incrementJetPt); + if(fIJMode) fIJHistosRecLeading->FillIntraJet( trackV, leadTrackV ); - fPhiCorrHistosAwayArea->FillTrackQA( partEta, - TVector2::Phi_mpi_pi( jet->Phi() - partPhi + TMath::Pi()), - partPt, - kTRUE); + delete trackV; } - } // end: phi-correlation - + + // ff and ij for background study + if(fBckgMode){ + if(fBckgType[0]!=-1) + FillBckgHistos(fBckgType[0], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, + fFFBckgHisto0RecCuts, fFFBckgHisto0RecLeading, + fIJBckgHisto0RecCuts, fIJBckgHisto0RecLeading, + fQABckgHisto0RecCuts); + if(fBckgType[1]!=-1) + FillBckgHistos(fBckgType[1], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, + fFFBckgHisto1RecCuts, fFFBckgHisto1RecLeading, + fIJBckgHisto1RecCuts, fIJBckgHisto1RecLeading, + fQABckgHisto1RecCuts); + if(fBckgType[2]!=-1) + FillBckgHistos(fBckgType[2], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, + fFFBckgHisto2RecCuts, fFFBckgHisto2RecLeading, + fIJBckgHisto2RecCuts, fIJBckgHisto2RecLeading, + fQABckgHisto2RecCuts); + if(fBckgType[3]!=-1) + FillBckgHistos(fBckgType[3], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, + fFFBckgHisto3RecCuts, fFFBckgHisto3RecLeading, + fIJBckgHisto3RecCuts, fIJBckgHisto3RecLeading, + fQABckgHisto3RecCuts); + if(fBckgType[4]!=-1) + FillBckgHistos(fBckgType[4], fTracksRecCuts, fJetsRecCuts, jet, leadTrackPt, leadTrackV, + fFFBckgHisto4RecCuts, fFFBckgHisto4RecLeading, + fIJBckgHisto4RecCuts, fIJBckgHisto4RecLeading, + fQABckgHisto4RecCuts); + } // end if(fBckgMode) + + + // phi correlation + if(fPhiCorrMode){ + for(Int_t it=0; itAt(it)); + + Float_t partEta = part->Eta(); + Float_t partPhi = part->Phi(); + Float_t partPt = part->Pt(); + + fPhiCorrHistosJetArea->FillTrackQA( partEta, + TVector2::Phi_mpi_pi( jet->Phi() - partPhi ), + partPt, + kTRUE); + + fPhiCorrHistosTransverseArea->FillTrackQA( partEta, + TVector2::Phi_mpi_pi( jet->Phi() - partPhi + TMath::Pi()/2), + partPt, + kTRUE); + + fPhiCorrHistosAwayArea->FillTrackQA( partEta, + TVector2::Phi_mpi_pi( jet->Phi() - partPhi + TMath::Pi()), + partPt, + kTRUE); + } + } // end: phi-correlation + + delete leadTrackV; + delete jettracklist; + + } // end: cut embedded ratio } // end: leading jet } // end: rec. jets after cuts @@ -4071,6 +4272,9 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) // secondaries if(fQAMode&1) FillSingleTrackHistosRecGen(0x0,fQATrackHistosSecRec,fTracksAODMCChargedSec,indexAODTrSec,isGenSec); + // high-pt occupancy effect + FillTwoTrackHistosRecGen(fTracksAODMCCharged, /*fTracksRecQualityCuts,indexAODTr,*/ isGenPrim); + // jet track eff Double_t sumPtGenLeadingJetRecEff = 0; @@ -4170,6 +4374,8 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) fJetsRecCuts->Clear(); fJetsGen->Clear(); fJetsRecEff->Clear(); + fJetsEmbedded->Clear(); + if(fBckgMode && (fBckgType[0]==kBckgClusters || fBckgType[1]==kBckgClusters || fBckgType[2]==kBckgClusters || fBckgType[3]==kBckgClusters || fBckgType[4]==kBckgClusters || @@ -4181,6 +4387,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) fBckgJetsGen->Clear(); } + //Post output data. PostData(1, fCommonHistList); } @@ -4273,15 +4480,41 @@ Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t t if(type==kTrackUndef) return 0; Int_t iCount = 0; - if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD){ - // all rec. tracks, esd filter mask, eta range + if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts || type==kTrackAODExtra || type==kTrackAODExtraonly){ + + TClonesArray *aodExtraTracks = dynamic_cast(fAOD->FindListObject("aodExtraTracks")); + if(!aodExtraTracks)return iCount; + for(int it =0; itGetEntries(); it++) { + AliVParticle *track = dynamic_cast ((*aodExtraTracks)[it]); + if (!track) continue; + + AliAODTrack *tr = dynamic_cast (track); + if(!tr)continue; + + if(type==kTrackAODExtraCuts || type==kTrackAODExtraonlyCuts){ + if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue; + + if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue; + if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue; + if(tr->Pt() < fTrackPtCut) continue; + } + + list->Add(tr); + iCount++; + } + } + + if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD || type==kTrackAODExtraCuts || type==kTrackAODExtra){ + + // all rec. tracks, esd filter mask, eta range for(Int_t it=0; itGetNumberOfTracks(); ++it){ AliAODTrack *tr = fAOD->GetTrack(it); - if(type == kTrackAODCuts || type==kTrackAODQualityCuts ){ + if(type == kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAODExtraCuts){ + if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask))) continue; if(type == kTrackAODCuts){ if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue; @@ -4490,7 +4723,7 @@ Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t typ if(!aodGenJets){ if(fDebug>0){ - if(fBranchGenJets.Length()) Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data()); + if(fBranchGenJets.Length()) Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data()); } if(fDebug>1)fAOD->Print(); return 0; @@ -4516,6 +4749,47 @@ Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t typ list->Sort(); return nGenJets; } + else if(type == kJetsEmbedded){ // embedded jets + + if(fBranchEmbeddedJets.Length()==0){ + Printf("%s:%d no embedded jet branch specified", (char*)__FILE__,__LINE__); + if(fDebug>1)fAOD->Print(); + return 0; + } + + TClonesArray *aodEmbeddedJets = 0; + if(fBranchEmbeddedJets.Length()) aodEmbeddedJets = dynamic_cast(fAODJets->FindListObject(fBranchEmbeddedJets.Data())); + if(!aodEmbeddedJets) aodEmbeddedJets = dynamic_cast(fAODJets->GetList()->FindObject(fBranchEmbeddedJets.Data())); + if(fAODExtension&&!aodEmbeddedJets) aodEmbeddedJets = dynamic_cast(fAODExtension->GetAOD()->FindListObject(fBranchEmbeddedJets.Data())); + + if(!aodEmbeddedJets){ + if(fBranchEmbeddedJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchEmbeddedJets.Data()); + if(fDebug>1)fAOD->Print(); + return 0; + } + + // Reorder jet pt and fill new temporary AliAODJet objects + Int_t nEmbeddedJets = 0; + + for(Int_t ij=0; ijGetEntries(); ++ij){ + + AliAODJet *tmp = dynamic_cast(aodEmbeddedJets->At(ij)); + if(!tmp) continue; + + if( tmp->Pt() < fJetPtCut ) continue; + if( tmp->Eta() < fJetEtaMin + || tmp->Eta() > fJetEtaMax + || tmp->Phi() < fJetPhiMin + || tmp->Phi() > fJetPhiMax ) continue; + + list->Add(tmp); + nEmbeddedJets++; + } + + list->Sort(); + + return nEmbeddedJets; + } else{ if(fDebug>0)Printf("%s:%d no such type %d",(char*)__FILE__,__LINE__,type); return 0; @@ -4767,6 +5041,7 @@ void AliAnalysisTaskFragmentationFunction::FillSingleTrackHistosRecGen(AliFragFu if(ptGen < fTrackPtCut) continue; if(trackQAGen) trackQAGen->FillTrackQA(etaGen, phiGen, ptGen); + Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track if(iRec>=0 && trackQARec) trackQARec->FillTrackQA(etaGen, phiGen, ptGen); @@ -4850,6 +5125,82 @@ void AliAnalysisTaskFragmentationFunction::FillJetTrackHistosRecGen(TObject* hi } } +// _____________________________________________________________________________________________________________________________________________ +void AliAnalysisTaskFragmentationFunction::FillTwoTrackHistosRecGen(TList* tracksGen, /*TList* tracksRec, const TArrayI& indexAODTr, */ const TArrayS& isRefGen){ + + Int_t nTracksGen = tracksGen->GetSize(); + + if(!nTracksGen) return; + + Int_t highPtIndices[nTracksGen]; + Int_t nHighPt = 0; + + for(Int_t iGen=0; iGen (tracksGen->At(iGen)); + if(!gentrack) continue; + Double_t ptGen = gentrack->Pt(); + Double_t etaGen = gentrack->Eta(); + Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi()); + + if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue; + if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue; + + if(ptGen>10 ){ + highPtIndices[nHighPt++] = iGen; + } + } + + + for(Int_t nHPT = 0; nHPT (tracksGen->At(indexHPT)); + if(!genHPTtrack) continue; + + Double_t trackMomHPT[3]; + genHPTtrack->PxPyPz(trackMomHPT); + TVector3 track3MomHPT(trackMomHPT); + + + Double_t distNN = 10; + Double_t ptNN = 0; + Int_t indexNN = -1; + + for(Int_t iGen=0; iGen (tracksGen->At(iGen)); + if(!gentrack) continue; + + Double_t ptGen = gentrack->Pt(); + Double_t etaGen = gentrack->Eta(); + Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi()); + + if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue; + if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue; + if(ptGen < fTrackPtCut) continue; + + + Double_t gentrackMom[3]; + gentrack->PxPyPz(gentrackMom); + TVector3 gentrack3Mom(gentrackMom); + + Double_t dR = gentrack3Mom.DeltaR(track3MomHPT); + + if(iGen != indexHPT && dRFill(entries); + + Double_t invPtGen = 0; + if(ptGen) invPtGen = 1/ptGen; + + Double_t invPtRec = 0; + if(ptRec) invPtRec = 1/ptRec; + + fh2SingleInvPtRecMnGenVsPtGen->Fill(ptGen,invPtRec - invPtGen); } } } @@ -4970,6 +5329,14 @@ void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxis(Float_t alp for (Int_t itrack=0; itrackGetSize(); itrack++){ + // embedded tracks + if( fUseExtraTracksBgr != 1){ + if(AliAODTrack* trackAOD = dynamic_cast (inputlist->At(itrack))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(inputlist->At(itrack)); if(!track)continue; Double_t trackMom[3]; @@ -5006,6 +5373,14 @@ void AliAnalysisTaskFragmentationFunction::GetTracksTiltedwrpJetAxisWindow(Float for (Int_t itrack=0; itrackGetSize(); itrack++) { + // embedded tracks + if( fUseExtraTracksBgr != 1){ + if(AliAODTrack* trackAOD = dynamic_cast (inputlist->At(itrack))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(inputlist->At(itrack)); if(!track)continue; Float_t trackEta = track->Eta(); @@ -5101,6 +5476,14 @@ void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJets(Int_t nCases, TLi // List of all tracks outside jet areas for (Int_t itrack=0; itrackGetSize(); itrack++){ + // embedded tracks + if( fUseExtraTracksBgr != 1){ + if(AliAODTrack* trackAOD = dynamic_cast (inputlist->At(itrack))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(inputlist->At(itrack)); if(!track)continue; Double_t trackMom[3]; @@ -5214,6 +5597,14 @@ void AliAnalysisTaskFragmentationFunction::GetTracksOutOfNJetsStat(Int_t nCases, for (Int_t itrack=0; itrackGetSize(); itrack++){ + // embedded tracks + if( fUseExtraTracksBgr != 1){ + if(AliAODTrack* trackAOD = dynamic_cast (inputlist->At(itrack))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(inputlist->At(itrack)); if(!track)continue; Double_t trackMom[3]; @@ -5299,6 +5690,14 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksOutOf1Jet(AliAODJet* for(Int_t it = 0; it (bgrCluster->GetTrack(it))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(bgrCluster->GetTrack(it)); if(!track) continue; @@ -5331,7 +5730,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range - if(!nBckgClusters) return; + if(nBckgClusters<3) return; // need at least 3 clusters (skipping 2 highest) Double_t* bgrDensity = new Double_t[nBckgClusters]; Int_t* indices = new Int_t[nBckgClusters]; @@ -5357,7 +5756,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl Double_t medianDensity = 0; if(TMath::Odd(nBckgClusters)){ - + Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters-1))]; medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex)); @@ -5378,7 +5777,7 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl Double_t clusterPt1 = medianCluster1->Pt(); Double_t area1 = medianCluster1->EffectiveAreaCharged(); if(area1>0) density1 = clusterPt1/area1; - + Double_t density2 = 0; Double_t clusterPt2 = medianCluster2->Pt(); Double_t area2 = medianCluster2->EffectiveAreaCharged(); @@ -5388,11 +5787,19 @@ void AliAnalysisTaskFragmentationFunction::GetClusterTracksMedian(TList* outputl medianCluster = ( (fRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 ); // select one randomly to avoid adding areas } - + Int_t nTracksJet = medianCluster->GetRefTracks()->GetEntries(); for(Int_t it = 0; it (medianCluster->GetTrack(it))){ + if(fUseExtraTracksBgr == 0 && ((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + if(fUseExtraTracksBgr == -1 && !((trackAOD->GetFlags() & AliESDtrack::kEmbedded)>0)) continue; + } + } + AliVParticle* track = dynamic_cast(medianCluster->GetTrack(it)); if(!track) continue; @@ -5450,13 +5857,13 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp if(type==kBckgOutLJ || type==kBckgOutAJ) { - TList* tracklistoutleading = new TList(); + TList* tracklistoutleading = new TList(); Double_t sumPtOutLeading = 0.; GetTracksOutOfNJets(1,inputtracklist, tracklistoutleading, inputjetlist, sumPtOutLeading); if(type==kBckgOutLJ) fh1OutLeadingMult->Fill(tracklistoutleading->GetSize()); for(Int_t it=0; itGetSize(); ++it){ - + AliVParticle* trackVP = (AliVParticle*)(tracklistoutleading->At(it)); if(!trackVP) continue; TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); @@ -5517,7 +5924,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp if(type==kBckgOutLJStat) fh1OutLeadingStatMult->Fill(tracklistoutleadingStat->GetSize()); for(Int_t it=0; itGetSize(); ++it){ - + AliVParticle* trackVP = dynamic_cast(tracklistoutleadingStat->At(it)); if(!trackVP) continue; TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); @@ -5621,7 +6028,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp fh1ASideMult->Fill(tracklistaside->GetSize()); for(Int_t it=0; itGetSize(); ++it){ - + AliVParticle* trackVP = (AliVParticle*)(tracklistaside->At(it)); if(!trackVP) continue; TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); @@ -5700,7 +6107,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp fh1PerpWindowMult->Fill(tracklistperpw->GetSize()); for(Int_t it=0; itGetSize(); ++it){ - + AliVParticle* trackVP = dynamic_cast(tracklistperpw->At(it)); if(!trackVP) continue; TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); @@ -5735,7 +6142,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp { if(type==kBckgOut2J) fh1Out2JetsMult->Fill(tracklistout2jets->GetSize()); for(Int_t it=0; itGetSize(); ++it){ - + AliVParticle* trackVP = dynamic_cast(tracklistout2jets->At(it)); if(!trackVP) continue; TLorentzVector* trackV = new TLorentzVector(trackVP->Px(),trackVP->Py(),trackVP->Pz(),trackVP->P()); @@ -5827,7 +6234,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp Float_t jetPt = jet->Pt(); Bool_t incrementJetPt = kTRUE; if(type==kBckgOut2JStat) - { + { if(fFFMode) ffbckghistocuts->FillFF( -1., jetPt, incrementJetPt, normFactor2Jets); if(fFFMode) ffbckghistoleading->FillFF( -1., leadTrackPt , incrementJetPt, normFactor2Jets); } @@ -5838,7 +6245,7 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp if(fFFMode) ffbckghistoleading->FillFF( -1., leadTrackPt , incrementJetPt, normFactor2Jets ); } } - + } if(type==kBckgOut3J || type==kBckgOutAJ) @@ -5894,8 +6301,6 @@ void AliAnalysisTaskFragmentationFunction::FillBckgHistos(Int_t type, TList* inp if(fFFMode) ffbckghistoleading->FillFF( -1., leadTrackPt , incrementJetPt ); } } - - } if(type==kBckgOut3JStat || type==kBckgOutAJStat) @@ -6050,7 +6455,4 @@ AliAODJet* AliAnalysisTaskFragmentationFunction::GetAODBckgSubJet(AliAODJet* jet } return tmpJet; - } - - diff --git a/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.h b/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.h index 7f528bf024c..36fe0019823 100644 --- a/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.h +++ b/PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.h @@ -340,12 +340,23 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { virtual void SetBranchRecBackClusters(const char* c){fBranchRecBckgClusters = c;} virtual void SetBranchGenJets(const char* c){fBranchGenJets = c;} virtual void SetBranchRecJets(const char* c){fBranchRecJets = c;} + virtual void SetBranchEmbeddedJets(const char* c){fBranchEmbeddedJets = c;} virtual void SetTrackCuts(Float_t trackPt = 0.15, Float_t trackEtaMin = -0.9, Float_t trackEtaMax = 0.9, Float_t trackPhiMin = 0., Float_t trackPhiMax = 2*TMath::Pi()) {fTrackPtCut = trackPt; fTrackEtaMin = trackEtaMin; fTrackEtaMax = trackEtaMax; fTrackPhiMin = trackPhiMin; fTrackPhiMax = trackPhiMax;} + virtual void UseExtraTracks() { fUseExtraTracks = 1;} + virtual void UseExtraonlyTracks() { fUseExtraTracks = -1;} + + virtual void UseExtraTracksBgr() { fUseExtraTracksBgr = 1;} + virtual void UseExtraonlyTracksBgr() { fUseExtraTracksBgr = -1;} + + virtual void SetCutFractionPtEmbedded(Float_t cut = 0) { fCutFractionPtEmbedded = cut; } + virtual void SetUseEmbeddedJetAxis(Bool_t b = kTRUE) { fUseEmbeddedJetAxis = b; } + virtual void SetUseEmbeddedJetPt(Bool_t b = kTRUE) { fUseEmbeddedJetPt = b; } + virtual void UseAODInputJets(Bool_t b) {fUseAODInputJets = b;} virtual void SetFilterMask(UInt_t i) {fFilterMask = i;} virtual void UsePhysicsSelection(Bool_t b) {fUsePhysicsSelection = b;} @@ -459,6 +470,8 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { void FillJetTrackHistosRecGen(TObject* histGen,TObject* histRec,Double_t jetPtGen,Double_t jetPtRec, TList* jetTrackList, const TList* tracksGen, const TArrayI& indexAODTr,const TArrayS& isRefGen, const Bool_t useRecJetPt); + void FillTwoTrackHistosRecGen(TList* tracksGen, /*TList* tracksRec, const TArrayI& indexAODTr, */ const TArrayS& isRefGen); + void FillSingleTrackResponse(THnSparse* hnResponse, TList* tracksGen, TList* tracksRec, const TArrayI& indexAODTr, const TArrayS& isGenPrim); void FillJetTrackResponse(THnSparse* hnResponsePt, THnSparse* hnResponseZ, THnSparse* hnResponseXi, Double_t jetPtGen, Double_t jetPtRec, TList* jetTrackList, @@ -476,10 +489,11 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { AliAODJet* GetAODBckgSubJet(AliAODJet* jet, Int_t method); // Consts - - enum {kTrackUndef=0, kTrackAOD, kTrackAODQualityCuts, kTrackAODCuts, kTrackKineAll, kTrackKineCharged, kTrackKineChargedAcceptance, + enum {kTrackUndef=0, kTrackAOD, kTrackAODQualityCuts, kTrackAODCuts, + kTrackAODExtra, kTrackAODExtraonly, kTrackAODExtraCuts, kTrackAODExtraonlyCuts, + kTrackKineAll, kTrackKineCharged, kTrackKineChargedAcceptance, kTrackAODMCAll, kTrackAODMCCharged, kTrackAODMCChargedAcceptance, kTrackAODMCChargedSec}; - enum {kJetsUndef=0, kJetsRec, kJetsRecAcceptance, kJetsGen, kJetsGenAcceptance, kJetsKine, kJetsKineAcceptance}; + enum {kJetsUndef=0, kJetsRec, kJetsRecAcceptance, kJetsGen, kJetsGenAcceptance, kJetsKine, kJetsKineAcceptance,kJetsEmbedded}; enum {kBckgPerp=0, kBckgOutLJ, kBckgOut2J, kBckgClusters, kBckgClustersOutLeading, kBckgOut3J, kBckgOutAJ, kBckgOutLJStat, kBckgOut2JStat, kBckgOut3JStat, kBckgOutAJStat, kBckgASide, kBckgASideWindow, kBckgPerpWindow}; @@ -502,7 +516,8 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { TString fBranchRecBackJets; // branch name for reconstructed background jets TString fBranchRecBckgClusters; // branch name for reconstructed background clusters TString fBranchGenJets; // branch name for generated jets - + TString fBranchEmbeddedJets; // branch name for embedded jets + Int_t fTrackTypeGen; // type of generated tracks Int_t fJetTypeGen; // type of generated jets @@ -521,6 +536,12 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { Float_t fTrackPhiMin; // track phi cut Float_t fTrackPhiMax; // track phi cut + Int_t fUseExtraTracks; // +/- 1: embedded extra/extra only tracks, default: 0 (ignore extra tracks) + Int_t fUseExtraTracksBgr; // +/- 1: background: use embedded extra/extra only tracks, default: 0 (ignore extra tracks) + Float_t fCutFractionPtEmbedded; // cut on ratio of embedded pt found in jet + Bool_t fUseEmbeddedJetAxis; // use axis of embedded jet for FF + Bool_t fUseEmbeddedJetPt; // use axis of embedded jet for FF + // jet cuts Float_t fJetPtCut; // jet transverse momentum cut Float_t fJetEtaMin; // jet eta cut @@ -564,6 +585,8 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { TList* fJetsRecCuts; //! jets from reonstructed tracks after jet cuts TList* fJetsGen; //! jets from generated tracks TList* fJetsRecEff; //! jets used for reconstruction efficiency histos + TList* fJetsEmbedded; //! jets used for embedding + TList* fBckgJetsRec; //! jets from reconstructed tracks TList* fBckgJetsRecCuts; //! jets from reonstructed tracks after jet cuts TList* fBckgJetsGen; //! jets from generated tracks @@ -753,11 +776,14 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { TH1F *fh1nRecJetsCuts; //! number of jets from reconstructed tracks per event TH1F *fh1nGenJets; //! number of jets from generated tracks per event TH1F *fh1nRecEffJets; //! number of jets for reconstruction eff per event + TH1F *fh1nEmbeddedJets; //! number of embedded jets per event + TH1F *fh1nRecBckgJetsCuts; //! number of jets from reconstructed tracks per event TH1F *fh1nGenBckgJets; //! number of jets from generated tracks per event TH2F *fh2PtRecVsGenPrim; //! association rec/gen MC: rec vs gen pt, primaries TH2F *fh2PtRecVsGenSec; //! association rec/gen MC: rec vs gen pt, secondaries + // tracking efficiency / secondaries AliFragFuncQATrackHistos* fQATrackHistosRecEffGen; //! tracking efficiency: generated primaries @@ -769,7 +795,9 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { AliFragFuncHistos* fFFHistosSecRec; //! secondary contamination: FF reconstructed secondaries // momentum resolution - THnSparse* fhnResponseSinglePt; //! single track response pt + THnSparse* fhnResponseSinglePt; //! single track response pt + TH2F* fh2SingleInvPtRecMnGenVsPtGen; //! single track response inv pt + THnSparse* fhnResponseJetTrackPt; //! jet track response pt THnSparse* fhnResponseJetZ; //! jet track response z THnSparse* fhnResponseJetXi; //! jet track response xi @@ -787,6 +815,19 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { TH1F *fh1MedianClustersMult; //! background multiplicity median cluster TH1F *fh1OutClustersMult; //! background multiplicity clusters outside leading jet + // embedding + TH1F* fh1FractionPtEmbedded; //! ratio embedded pt in rec jet to embedded jet pt + TH1F* fh1IndexEmbedded; //! index embedded jet matching to leading rec jet + TH2F* fh2DeltaPtVsJetPtEmbedded; //! delta pt rec - embedded jet + TH2F* fh2DeltaPtVsRecJetPtEmbedded; //! delta pt rec - embedded jet + TH1F* fh1DeltaREmbedded; //! delta R rec - embedded jet + + // 2track resolution + TH2F* fh2ptVsDistNN_pt50_rec; //! + TH2F* fh2ptVsDistNN_pt50_nonRec; //! + TH2F* fh2ptVsDistNN_pt10_rec; //! + TH2F* fh2ptVsDistNN_pt10_nonRec; //! + AliFragFuncQATrackHistos* fQABckgHisto0RecCuts; //! track QA: reconstructed tracks after cuts AliFragFuncQATrackHistos* fQABckgHisto0Gen; //! track QA: generated tracks AliFragFuncQATrackHistos* fQABckgHisto1RecCuts; //! track QA: reconstructed tracks after cuts @@ -844,7 +885,7 @@ class AliAnalysisTaskFragmentationFunction : public AliAnalysisTaskSE { TRandom3* fRandom; // TRandom3 for background estimation Int_t fBckgSubMethod; // Bckg method: 1 = leading jet excluded, 2 = 2 most energetic jets excluded - ClassDef(AliAnalysisTaskFragmentationFunction, 10); + ClassDef(AliAnalysisTaskFragmentationFunction, 11); }; #endif diff --git a/PWG4/JetTasks/AliAnalysisTaskJetChem.cxx b/PWG4/JetTasks/AliAnalysisTaskJetChem.cxx index 7161b2f4e48..ffc4f9c2e47 100644 --- a/PWG4/JetTasks/AliAnalysisTaskJetChem.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskJetChem.cxx @@ -228,12 +228,7 @@ AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem() { // destructor - - if(fTracksRecCuts) delete fTracksRecCuts; - if(fJetsRecCuts) delete fJetsRecCuts; if(fListK0s) delete fListK0s; - if(fCommonHistList) delete fCommonHistList; - } //________________________________________________________________________________________________________________________________ diff --git a/PWG4/JetTasks/AliFragmentationFunctionCorrections.cxx b/PWG4/JetTasks/AliFragmentationFunctionCorrections.cxx new file mode 100644 index 00000000000..c061320d501 --- /dev/null +++ b/PWG4/JetTasks/AliFragmentationFunctionCorrections.cxx @@ -0,0 +1,4013 @@ +// ************************************************************************* +// * * +// * corrections to Fragmentation Functions * +// * * +// ************************************************************************* + + +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id: */ + +#include "TMath.h" +#include "TH2F.h" +#include "THnSparse.h" +#include "TFile.h" +#include "TDirectory.h" +#include "AliCFUnfolding.h" +#include "AliFragmentationFunctionCorrections.h" + +/////////////////////////////////////////////////////////////////////////////// +// // +// correction class for ouput of AliAnalysisTaskFragmentationFunction // +// // +// corrections for: reconstruction efficiency, momentum resolution, // +// secondaries, UE / HI background, fluctuations // +// back-correction on jet energy on dN/dxi // +// // +// read MC ouput and write out efficiency histos, response matrices etc. // +// read measured distributions and apply efficiency, response matrices etc. // +// // +// contact: o.busch@gsi.de // +// // +/////////////////////////////////////////////////////////////////////////////// + + +ClassImp(AliFragmentationFunctionCorrections) + +//________________________________________________________________________ +AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections() + : TObject() + ,fDebug(0) + ,fNJetPtSlices(0) + ,fJetPtSlices(0) + ,fNJets(0) + ,fNJetsBgr(0) + ,fNHistoBinsSinglePt(0) + ,fHistoBinsSinglePt(0) + ,fNHistoBinsPt(0) + ,fNHistoBinsZ(0) + ,fNHistoBinsXi(0) + ,fHistoBinsPt(0) + ,fHistoBinsZ(0) + ,fHistoBinsXi(0) + ,fNCorrectionLevels(0) + ,fCorrFF(0) + ,fNCorrectionLevelsBgr(0) + ,fCorrBgr(0) + ,fNCorrectionLevelsSinglePt(0) + ,fCorrSinglePt(0) + ,fh1FFXiShift(0) + ,fh1EffSinglePt(0) + ,fh1EffPt(0) + ,fh1EffZ(0) + ,fh1EffXi(0) + ,fh1EffBgrPt(0) + ,fh1EffBgrZ(0) + ,fh1EffBgrXi(0) + ,fh1FFTrackPtBackFolded(0) + ,fh1FFZBackFolded(0) + ,fh1FFXiBackFolded(0) + ,fh1FFRatioTrackPtFolded(0) + ,fh1FFRatioZFolded(0) + ,fh1FFRatioXiFolded(0) + ,fh1FFRatioTrackPtBackFolded(0) + ,fh1FFRatioZBackFolded(0) + ,fh1FFRatioXiBackFolded(0) + ,fh1SingleTrackPtBackFolded(0) + ,fh1RatioSingleTrackPtFolded(0) + ,fh1RatioSingleTrackPtBackFolded(0) + ,fhnResponseSinglePt(0) + ,fhnResponsePt(0) + ,fhnResponseZ(0) + ,fhnResponseXi(0) + ,fh1FFTrackPtPrior(0) + ,fh1FFZPrior(0) + ,fh1FFXiPrior(0) + ,fh1SecCorrSinglePt(0) +{ + // default constructor +} + +//________________________________________________________________________________________________________________________ +AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections ©) + : TObject() + ,fDebug(copy.fDebug) + ,fNJetPtSlices(copy.fNJetPtSlices) + ,fJetPtSlices(copy.fJetPtSlices) + ,fNJets(copy.fNJets) + ,fNJetsBgr(copy.fNJetsBgr) + ,fNHistoBinsSinglePt(fNHistoBinsSinglePt) + ,fHistoBinsSinglePt(fHistoBinsSinglePt) + ,fNHistoBinsPt(copy.fNHistoBinsPt) + ,fNHistoBinsZ(copy.fNHistoBinsZ) + ,fNHistoBinsXi(copy.fNHistoBinsXi) + ,fHistoBinsPt(copy.fHistoBinsPt) + ,fHistoBinsZ(copy.fHistoBinsZ) + ,fHistoBinsXi(copy.fHistoBinsXi) + ,fNCorrectionLevels(copy.fNCorrectionLevels) + ,fCorrFF(copy.fCorrFF) + ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr) + ,fCorrBgr(copy.fCorrBgr) + ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt) + ,fCorrSinglePt(copy.fCorrSinglePt) + ,fh1FFXiShift(copy.fh1FFXiShift) + ,fh1EffSinglePt(fh1EffSinglePt) + ,fh1EffPt(copy.fh1EffPt) + ,fh1EffZ(copy.fh1EffZ) + ,fh1EffXi(copy.fh1EffXi) + ,fh1EffBgrPt(copy.fh1EffBgrPt) + ,fh1EffBgrZ(copy.fh1EffBgrZ) + ,fh1EffBgrXi(copy.fh1EffBgrXi) + ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded) + ,fh1FFZBackFolded(copy.fh1FFZBackFolded) + ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded) + ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded) + ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded) + ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded) + ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded) + ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded) + ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded) + ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded) + ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded) + ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded) + ,fhnResponseSinglePt(copy.fhnResponseSinglePt) + ,fhnResponsePt(copy.fhnResponsePt) + ,fhnResponseZ(copy.fhnResponseZ) + ,fhnResponseXi(copy.fhnResponseXi) + ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior) + ,fh1FFZPrior(copy.fh1FFZPrior) + ,fh1FFXiPrior(copy.fh1FFXiPrior) + ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt) +{ + // copy constructor + +} + +// ______________________________________________________________________________________________________________________________ +AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o) +{ + // assignment + + if(this!=&o){ + TObject::operator=(o); + fDebug = o.fDebug; + fNJetPtSlices = o.fNJetPtSlices; + fJetPtSlices = o.fJetPtSlices; + fNJets = o.fNJets; + fNJetsBgr = o.fNJetsBgr; + fNHistoBinsSinglePt = o.fNHistoBinsSinglePt; + fHistoBinsSinglePt = o.fHistoBinsSinglePt; + fNHistoBinsPt = o.fNHistoBinsPt; + fNHistoBinsZ = o.fNHistoBinsZ; + fNHistoBinsXi = o.fNHistoBinsXi; + fHistoBinsPt = o.fHistoBinsPt; + fHistoBinsZ = o.fHistoBinsZ; + fHistoBinsXi = o.fHistoBinsXi; + fNCorrectionLevels = o.fNCorrectionLevels; + fCorrFF = o.fCorrFF; + fNCorrectionLevelsBgr = o.fNCorrectionLevelsBgr; + fCorrBgr = o.fCorrBgr; + fNCorrectionLevelsSinglePt = o.fNCorrectionLevelsSinglePt; + fCorrSinglePt = o.fCorrSinglePt; + fh1FFXiShift = o.fh1FFXiShift; + fh1EffSinglePt = o.fh1EffSinglePt; + fh1EffPt = o.fh1EffPt; + fh1EffZ = o.fh1EffZ; + fh1EffXi = o.fh1EffXi; + fh1EffBgrPt = o.fh1EffBgrPt; + fh1EffBgrZ = o.fh1EffBgrZ; + fh1EffBgrXi = o.fh1EffBgrXi; + fh1FFTrackPtBackFolded = o.fh1FFTrackPtBackFolded; + fh1FFZBackFolded = o.fh1FFZBackFolded; + fh1FFXiBackFolded = o.fh1FFXiBackFolded; + fh1FFRatioTrackPtFolded = o.fh1FFRatioTrackPtFolded; + fh1FFRatioZFolded = o.fh1FFRatioZFolded; + fh1FFRatioXiFolded = o.fh1FFRatioXiFolded; + fh1FFRatioTrackPtBackFolded = o.fh1FFRatioTrackPtBackFolded; + fh1FFRatioZBackFolded = o.fh1FFRatioZBackFolded; + fh1FFRatioXiBackFolded = o.fh1FFRatioXiBackFolded; + fh1SingleTrackPtBackFolded = o.fh1SingleTrackPtBackFolded; + fh1RatioSingleTrackPtFolded = o.fh1RatioSingleTrackPtFolded; + fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded; + fhnResponseSinglePt = o.fhnResponseSinglePt; + fhnResponsePt = o.fhnResponsePt; + fhnResponseZ = o.fhnResponseZ; + fhnResponseXi = o.fhnResponseXi; + fh1FFTrackPtPrior = o.fh1FFTrackPtPrior; + fh1FFZPrior = o.fh1FFZPrior; + fh1FFXiPrior = o.fh1FFXiPrior; + fh1SecCorrSinglePt = o.fh1SecCorrSinglePt; + } + + return *this; +} + +//_________________________________________________________________________ +AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections() +{ + // destructor + + if(fJetPtSlices) delete fJetPtSlices; + if(fNJets) delete fNJets; + if(fNJetsBgr) delete fNJetsBgr; + + DeleteHistoArray(fh1FFXiShift); + + DeleteHistoArray(fh1EffPt); + DeleteHistoArray(fh1EffXi); + DeleteHistoArray(fh1EffZ ); + + DeleteHistoArray(fh1EffBgrPt); + DeleteHistoArray(fh1EffBgrXi); + DeleteHistoArray(fh1EffBgrZ); + + // unfolding + + DeleteHistoArray(fh1FFTrackPtBackFolded); + DeleteHistoArray(fh1FFZBackFolded); + DeleteHistoArray(fh1FFXiBackFolded); + + DeleteHistoArray(fh1FFRatioTrackPtFolded); + DeleteHistoArray(fh1FFRatioZFolded); + DeleteHistoArray(fh1FFRatioXiFolded); + + DeleteHistoArray(fh1FFRatioTrackPtBackFolded); + DeleteHistoArray(fh1FFRatioZBackFolded); + DeleteHistoArray(fh1FFRatioXiBackFolded); + + DeleteTHnSparseArray(fhnResponsePt); + DeleteTHnSparseArray(fhnResponseZ); + DeleteTHnSparseArray(fhnResponseXi); + + DeleteHistoArray(fh1FFTrackPtPrior); + DeleteHistoArray(fh1FFZPrior); + DeleteHistoArray(fh1FFXiPrior); + + + // clean up corrected FF + + for(Int_t c=0; c=fArraySize){ + Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__); + return; + } + + if(histPt){ + TString name = histPt->GetName(); + if(fCorrLabel.Length()>0) name += "_"+fCorrLabel; + histPt->SetName(name); + + if(!(histPt->GetSumw2())) histPt->Sumw2(); + + new(fh1CorrFFTrackPt[slice]) TH1F(*histPt); + } + + if(histZ){ + TString name = histZ->GetName(); + if(fCorrLabel.Length()>0) name += "_"+fCorrLabel; + histZ->SetName(name); + + if(!(histZ->GetSumw2())) histZ->Sumw2(); + + new(fh1CorrFFZ[slice]) TH1F(*histZ); + } + + if(histXi){ + TString name = histXi->GetName(); + if(fCorrLabel.Length()>0) name += "_"+fCorrLabel; + histXi->SetName(name); + + if(!(histXi->GetSumw2())) histXi->Sumw2(); + + new(fh1CorrFFXi[slice]) TH1F(*histXi); + } +} + +//___________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi) +{ + // replace histo in array + + if(slice>=fArraySize){ + Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__); + return; + } + + if(histPt){ + if(!(histPt->GetSumw2())) histPt->Sumw2(); + + TString name = histPt->GetName(); + histPt->SetName(name); + + delete fh1CorrFFTrackPt[slice]; + fh1CorrFFTrackPt[slice] = new TH1F; + new(fh1CorrFFTrackPt[slice]) TH1F(*histPt); + } + + if(histZ){ + if(!(histZ->GetSumw2())) histZ->Sumw2(); + + TString name = histZ->GetName(); + histZ->SetName(name); + + delete fh1CorrFFZ[slice]; + fh1CorrFFZ[slice] = new TH1F; + new(fh1CorrFFZ[slice]) TH1F(*histZ); + } + + if(histXi){ + if(!(histXi->GetSumw2())) histXi->Sumw2(); + + TString name = histXi->GetName(); + histXi->SetName(name); + + delete fh1CorrFFXi[slice]; + fh1CorrFFXi[slice] = new TH1F; + new(fh1CorrFFXi[slice]) TH1F(*histXi); + } +} + +// ___________________________________________________________________________________________ +TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice) +{ + // return pt histo + + if(slice>=fArraySize){ + Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__); + return 0; + } + + return fh1CorrFFTrackPt[slice]; + +} + +// ______________________________________________________________________________________ +TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice) +{ + // return z histo + + if(slice>=fArraySize){ + Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__); + return 0; + } + + return fh1CorrFFZ[slice]; +} + +// ________________________________________________________________________________________ +TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice) +{ + // return xi histo + + if(slice>=fArraySize){ + Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__); + return 0; + } + + return fh1CorrFFXi[slice]; +} + +// __________________________________________________________________________ +void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const +{ + // delete array of TH1 + + if(hist){ + for(Int_t i=0; i= fgMaxNCorrectionLevels){ + Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__); + return; + } + + fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices); + fNCorrectionLevels++; +} + + +//________________________________________________________________________________ +void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label) +{ + // increase corr level for bgr FF + + if(!fNJetPtSlices){ + if(fDebug>0) Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__); + return; + } + + if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){ + Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__); + return; + } + + fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices); + fNCorrectionLevelsBgr++; +} + +//_____________________________________________________________________________________ +void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label) +{ + // increase corr level single pt spec + + Int_t nJetPtSlicesSingle = 1; // pro forma + + if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){ + Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__); + return; + } + + fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle); + fNCorrectionLevelsSinglePt++; +} + +//_____________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices) +{ + // set jet pt slices + // once slices are known, can book all other histos + + fNJetPtSlices = nJetPtSlices; + + const Int_t size = nJetPtSlices+1; + fJetPtSlices = new TArrayF(size,bins); + + // nJets array + + fNJets = new TArrayF(size); + fNJets->Reset(0); + + fNJetsBgr = new TArrayF(size); + fNJetsBgr->Reset(0); + + // histos + + fNCorrectionLevels = 0; + fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels]; + AddCorrectionLevel(); // first 'correction' level = raw FF + + fNCorrectionLevelsBgr = 0; + fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels]; + AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist + + fh1FFXiShift = BookHistoArray(); + + // eff histos + + fh1EffPt = BookHistoArray(); + fh1EffXi = BookHistoArray(); + fh1EffZ = BookHistoArray(); + + fh1EffBgrPt = BookHistoArray(); + fh1EffBgrXi = BookHistoArray(); + fh1EffBgrZ = BookHistoArray(); + + + // unfolding + + fh1FFTrackPtBackFolded = BookHistoArray(); + fh1FFXiBackFolded = BookHistoArray(); + fh1FFZBackFolded = BookHistoArray(); + + fh1FFRatioTrackPtFolded = BookHistoArray(); + fh1FFRatioXiFolded = BookHistoArray(); + fh1FFRatioZFolded = BookHistoArray(); + + fh1FFRatioTrackPtBackFolded = BookHistoArray(); + fh1FFRatioXiBackFolded = BookHistoArray(); + fh1FFRatioZBackFolded = BookHistoArray(); + + // + + fhnResponsePt = BookTHnSparseArray(); + fhnResponseZ = BookTHnSparseArray(); + fhnResponseXi = BookTHnSparseArray(); + + fh1FFTrackPtPrior = BookHistoArray(); + fh1FFZPrior = BookHistoArray(); + fh1FFXiPrior = BookHistoArray(); + + + // histos bins + + fNHistoBinsPt = new Int_t[fNJetPtSlices]; + fNHistoBinsXi = new Int_t[fNJetPtSlices]; + fNHistoBinsZ = new Int_t[fNJetPtSlices]; + + for(Int_t i=0; i=fNJetPtSlices){ + Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice); + return; + } + + if(type == kFlagPt){ + fNHistoBinsPt[jetPtSlice] = sizeBins-1; + fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins); + } + else if(type==kFlagZ){ + fNHistoBinsZ[jetPtSlice] = sizeBins-1; + fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins); + } + + else if(type==kFlagXi){ + fNHistoBinsXi[jetPtSlice] = sizeBins-1; + fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins); + } +} + +//__________________________________________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type) +{ + // set histo bins for jet pt slice + // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...) + // array size of binsLimits: nBinsLimits + // array size of binsWidth: nBinsLimits-1 + // binsLimits have to be in increasing order + // if binning undefined for any slice, original binning will be used + + if(!fNJetPtSlices){ + Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__); + return; + } + + if(jetPtSlice>=fNJetPtSlices){ + Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice); + return; + } + + + Double_t binLimitMin = binsLimits[0]; + Double_t binLimitMax = binsLimits[nBinsLimits-1]; + + Double_t binLimit = binLimitMin; // start value + + Int_t sizeUpperLim = 10000; //static_cast(binLimitMax/binsWidth[0])+1; + TArrayD binsArray(sizeUpperLim); + Int_t nBins = 0; + binsArray.SetAt(binLimitMin,nBins++); + + while(binLimit= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor + } + + Double_t currentBinWidth = binsWidth[currentSlice]; + binLimit += currentBinWidth; + + binsArray.SetAt(binLimit,nBins++); + } + + Double_t* bins = binsArray.GetArray(); + + SetHistoBins(jetPtSlice,nBins,bins,type); +} + +//__________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins) +{ + // set histo bins for inclusive pt spectra + // if binning undefined, original binning will be used + + fNHistoBinsSinglePt = sizeBins-1; + + fHistoBinsSinglePt = new TArrayD(sizeBins); + fHistoBinsSinglePt->Set(sizeBins,bins); +} + +//__________________________________________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth) +{ + // set histo bins for inclusive pt spectra + // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...) + // array size of binsLimits: nBinsLimits + // array size of binsWidth: nBinsLimits-1 + // binsLimits have to be in increasing order + // if binning undefined for any slice, original binning will be used + + + Double_t binLimitMin = binsLimits[0]; + Double_t binLimitMax = binsLimits[nBinsLimits-1]; + + Double_t binLimit = binLimitMin; // start value + + Int_t sizeUpperLim = 10000; //static_cast(binLimitMax/binsWidth[0])+1; + TArrayD binsArray(sizeUpperLim); + Int_t nBins = 0; + binsArray.SetAt(binLimitMin,nBins++); + + while(binLimit= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor + } + + Double_t currentBinWidth = binsWidth[currentSlice]; + binLimit += currentBinWidth; + + binsArray.SetAt(binLimit,nBins++); + } + + Double_t* bins = binsArray.GetArray(); + + SetHistoBinsSinglePt(nBins,bins); +} + +//____________________________________________________________________________________ +void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets) +{ + // FF normalization: divide by bin width and normalize to entries/jets + // should also work for TH2, but to be tested !!! + + if(nJets == 0){ // nothing to do + if(fDebug>0) Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__); + return; + } + + Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ(); + + for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?) + + Double_t binWidth = hist->GetBinWidth(bin); + Double_t binCont = hist->GetBinContent(bin); + Double_t binErr = hist->GetBinError(bin); + + binCont /= binWidth; + binErr /= binWidth; + + hist->SetBinContent(bin,binCont); + hist->SetBinError(bin,binErr); + } + + hist->Scale(1/nJets); +} + +//_____________________________________________________ +void AliFragmentationFunctionCorrections::NormalizeFF() +{ + // normalize FF + + if(fNCorrectionLevels>1){ + Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__); + } + + for(Int_t i=0; i0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i)); + + NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF + NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i)); + NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i)); + } +} + +//______________________________________________________ +void AliFragmentationFunctionCorrections::NormalizeBgr() +{ + // normalize bgr/UE FF + + if(fNCorrectionLevelsBgr>1){ + Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__); + } + + for(Int_t i=0; iGetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF + NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i)); + NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i)); + } + +} + +//__________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID) +{ + // read raw FF - standard dir/list name + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + ReadRawFF(strfile,strdir,strlist,strFFID); +} + +//____________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID) +{ + // get raw FF from input file, project in jet pt slice + // normalization done separately + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read FF from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + gDirectory->cd(strdir); + + TList* list = 0; + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + + TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); + TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data())); + TString hnameZ(Form("fh2FFZ%s",strFFID.Data())); + TString hnameXi(Form("fh2FFXi%s",strFFID.Data())); + + TH1F* fh1FFJetPt = (TH1F*) list->FindObject(hnameJetPt); + TH2F* fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt); + TH2F* fh2FFZ = (TH2F*) list->FindObject(hnameZ); + TH2F* fh2FFXi = (TH2F*) list->FindObject(hnameXi); + + if(!fh1FFJetPt) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; } + if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; } + if(!fh2FFZ) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data()); return; } + if(!fh2FFXi) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data()); return; } + + fh1FFJetPt->SetDirectory(0); + fh2FFTrackPt->SetDirectory(0); + fh2FFZ->SetDirectory(0); + fh2FFXi->SetDirectory(0); + + f.Close(); + + + // nJets per bin + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh1FFJetPt->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh1FFJetPt->FindBin(jetPtUpLim)) - 1; + + Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp); + + fNJets->SetAt(nJetsBin,i); + + if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast(jetPtLoLim),static_cast(jetPtUpLim),fNJets->At(i)); + } + + // projections: FF + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1; + + if(binUp > fh2FFTrackPt->GetNbinsX()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + // appendix 'unbinned' to avoid histos with same name after rebinning + TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + TH1F* projZ = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o"); + TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o"); + + if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray()); + + projPt->SetNameTitle(strNameFFPt,""); + projZ->SetNameTitle(strNameFFZ,""); + projXi->SetNameTitle(strNameFFXi,""); + + // raw FF = corr level 0 + fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi); + } +} + +//_____________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID) +{ + // read raw FF - standard dir/list name + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID); +} + +//_______________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID) +{ + // get raw FF from input file, project in jet pt slice + // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID + // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts" + // normalization done separately + + TString strID = strBgrID + strFFID; + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data()); + + gDirectory->cd(strdir); + + TList* list = 0; + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + + TString hnameNJets = "fh1nRecJetsCuts"; + TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization + TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data())); + TString hnameBgrZ(Form("fh2FFZ%s",strID.Data())); + TString hnameBgrXi(Form("fh2FFXi%s",strID.Data())); + + TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets + TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt); + TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt); + TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ); + TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi); + + if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; } + if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; } + if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; } + if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; } + if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; } + + fh1FFJetPtBgr->SetDirectory(0); + fh1NJets->SetDirectory(0); + fh2FFTrackPtBgr->SetDirectory(0); + fh2FFZBgr->SetDirectory(0); + fh2FFXiBgr->SetDirectory(0); + + f.Close(); + + // nJets per bin + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh1FFJetPtBgr->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1; + + Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp); + Double_t scaleF = 1; + + //if(strBgrID.Contains("Out2Jets")){ // scale by ratio 2 jets events / all events + // scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX()) + // / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());} + + + if(strBgrID.Contains("OutAllJets")){ // scale by ratio >3 jets events / all events + scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX()) + / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX()); + } + + fNJetsBgr->SetAt(nJetsBin*scaleF,i); + + if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f", + static_cast(jetPtLoLim),static_cast(jetPtUpLim),nJetsBin,scaleF); + + } + + // projections: FF + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1; + + if(binUp > fh2FFTrackPtBgr->GetNbinsX()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + // appendix 'unbinned' to avoid histos with same name after rebinning + TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o"); + TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o"); + + if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray()); + + projPt->SetNameTitle(strNameBgrPt,""); + projZ->SetNameTitle(strNameBgrZ,""); + projXi->SetNameTitle(strNameBgrXi,""); + + // raw bgr = corr level 0 + fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi); + } +} + +//_____________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID) +{ + // read raw FF - standard dir/list name + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID); +} + +//_______________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID) +{ + // get raw FF from input file, project in jet pt slice + // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts" + // normalization done separately + + TString strBgrID = "BckgEmbed"; + TString strID = strBgrID + strFFID; + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data()); + + gDirectory->cd(strdir); + + TList* list = 0; + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + + TString hnameNJets = "fh1nRecJetsCuts"; + TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); + TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data())); + TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data())); + TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data())); + + TH1F* fh1NJets = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets + TH1F* fh1FFJetPtBgr = (TH1F*) list->FindObject(hnameJetPt); + TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt); + TH2F* fh2FFZBgr = (TH2F*) list->FindObject(hnameBgrZ); + TH2F* fh2FFXiBgr = (TH2F*) list->FindObject(hnameBgrXi); + + if(!fh1FFJetPtBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data()); return; } + if(!fh1NJets) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data()); return; } + if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; } + if(!fh2FFZBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data()); return; } + if(!fh2FFXiBgr) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data()); return; } + + fh1FFJetPtBgr->SetDirectory(0); + fh1NJets->SetDirectory(0); + fh2FFTrackPtBgr->SetDirectory(0); + fh2FFZBgr->SetDirectory(0); + fh2FFXiBgr->SetDirectory(0); + + f.Close(); + + // nJets per bin + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh1FFJetPtBgr->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1; + + Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp); + Double_t scaleF = 1; + + fNJetsBgr->SetAt(nJetsBin*scaleF,i); + + if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f", + static_cast(jetPtLoLim),static_cast(jetPtUpLim),nJetsBin,scaleF); + + } + + // projections: FF + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1; + + if(binUp > fh2FFTrackPtBgr->GetNbinsX()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + // appendix 'unbinned' to avoid histos with same name after rebinning + TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + TH1F* projZ = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o"); + TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o"); + + if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) projZ = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray()); + + projPt->SetNameTitle(strNameBgrPt,""); + projZ->SetNameTitle(strNameBgrZ,""); + projXi->SetNameTitle(strNameBgrXi,""); + + // raw bgr = corr level 0 + fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi); + } +} + + +//__________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile) +{ + // write histos to file + // skip histos with 0 entries + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile f(strfile,outfileOption); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + if(strdir && strdir.Length()){ + TDirectory* dir = f.mkdir(strdir); + dir->cd(); + } + + for(Int_t i=0; iGetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write(); + for(Int_t c=0; cGetZ(i)->GetEntries()) fCorrFF[c]->GetZ(i)->Write(); + for(Int_t c=0; cGetXi(i)->GetEntries()) fCorrFF[c]->GetXi(i)->Write(); + + if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write(); + + for(Int_t c=0; cGetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write(); + for(Int_t c=0; cGetZ(i)->GetEntries()) fCorrBgr[c]->GetZ(i)->Write(); + for(Int_t c=0; cGetXi(i)->GetEntries()) fCorrBgr[c]->GetXi(i)->Write(); + + + if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write(); + if(fh1FFZBackFolded[i] && fh1FFZBackFolded[i]->GetEntries()) fh1FFZBackFolded[i]->Write(); + if(fh1FFXiBackFolded[i] && fh1FFXiBackFolded[i]->GetEntries()) fh1FFXiBackFolded[i]->Write(); + + + if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write(); + if(fh1FFRatioZFolded[i] && fh1FFRatioZFolded[i]->GetEntries()) fh1FFRatioZFolded[i]->Write(); + if(fh1FFRatioXiFolded[i] && fh1FFRatioXiFolded[i]->GetEntries()) fh1FFRatioXiFolded[i]->Write(); + + if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write(); + if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries()) fh1FFRatioZBackFolded[i]->Write(); + if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries()) fh1FFRatioXiBackFolded[i]->Write(); + + } + + // inclusive track pt + + for(Int_t c=0; cGetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write(); + if(fh1SingleTrackPtBackFolded) fh1SingleTrackPtBackFolded->Write(); + if(fh1RatioSingleTrackPtFolded) fh1RatioSingleTrackPtFolded->Write(); + if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write(); + + f.Close(); +} + +//____________________________________________________________________________________________________________________________________ +THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst) +{ + // copy 1-dimensional histo to THnSparse + // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1 + // histos with variable bin size are supported + + // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user + + THnSparseF* fhnHist; + + Int_t nBins = hist->GetXaxis()->GetNbins(); + Int_t nBinsVec[1] = { nBins }; + + const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray(); + + if(binsVec){ // binsVec only neq NULL if histo was rebinned before + + fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/); + fhnHist->SetBinEdges(0,binsVec); + } + else{ // uniform bin size + + Double_t xMin = hist->GetXaxis()->GetXmin(); + Double_t xMax = hist->GetXaxis()->GetXmax(); + + Double_t binMinVec[1] = { xMin }; + Double_t binMaxVec[1] = { xMax }; + + fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec); + } + + + for(Int_t bin=1; bin<=nBins; bin++){ + + Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin); + + Double_t binCoord[] = {binCenter}; + fhnHist->Fill(binCoord,1); // initially need to create the bin + + Long64_t binIndex = fhnHist->GetBin(binCoord); + + Double_t cont = hist->GetBinContent(bin); + Double_t err = hist->GetBinError(bin); + + if(fillConst){ + cont = 1; + err = 0; + } + + fhnHist->SetBinContent(binIndex,cont); + fhnHist->SetBinError(binIndex,err); + } + + return fhnHist; +} + +//______________________________________________________________________________________________________________________________________________ +THnSparse* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter, + const Bool_t useCorrelatedErrors, const THnSparse* hnPrior) +{ + // unfold input histo + + AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior); // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0) + unfolding.SetMaxNumberOfIterations(nIter); + // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!! + // if(useSmoothing) unfolding.UseSmoothing(); + + if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors(); + unfolding.Unfold(); + + THnSparse* unfolded = unfolding.GetUnfolded(); + + TString hnameUnf = hnHist->GetName(); + hnameUnf.Append("_unf"); + unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle()); + + return unfolded; +} + +//___________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type) +{ + // loop over jet pt slices and unfold dN/dpt spectra + + TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel(); + if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold"); + + for(Int_t i=0; iGetTrackPt(i); // level -2: before unfolding, level -1: unfolded + if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i); // level -2: before unfolding, level -1: unfolded + if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i); // level -2: before unfolding, level -1: unfolded + + THnSparse* hnResponse = 0; + if(type == kFlagPt) hnResponse = fhnResponsePt[i]; + if(type == kFlagZ) hnResponse = fhnResponseZ[i]; + if(type == kFlagXi) hnResponse = fhnResponseXi[i]; + + TH1F* hPrior = 0; + if(type == kFlagPt && fh1FFTrackPtPrior[i] && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i]; + if(type == kFlagZ && fh1FFZPrior[i] && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFZPrior[i]; + if(type == kFlagXi && fh1FFXiPrior[i] && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFXiPrior[i]; + + TString histNameTHn = hist->GetName(); + histNameTHn.ReplaceAll("TH1","THn"); + + TString priorNameTHn; + if(hPrior){ + priorNameTHn = hPrior->GetName(); + priorNameTHn.ReplaceAll("TH1","THn"); + } + + TString histNameBackFolded = hist->GetName(); + histNameBackFolded.Append("_backfold"); + + TString histNameRatioFolded = hist->GetName(); + histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF"); + histNameRatioFolded.Append("_unfold"); + + TString histNameRatioBackFolded = hist->GetName(); + histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF"); + histNameRatioBackFolded.Append("_backfold"); + + THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle()); + THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff + THnSparse* hnPrior = 0; + if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle()); + + THnSparse* hnUnfolded + = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior); + + TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); + hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle()); + + if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0); + if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0); + if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded); + + // backfolding: apply response matrix to unfolded spectrum + TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); + hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle()); + + if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded; + if(type == kFlagZ) fh1FFZBackFolded[i] = hBackFolded; + if(type == kFlagXi) fh1FFXiBackFolded[i] = hBackFolded; + + // ratio unfolded to original histo + TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded); + hRatioUnfolded->Reset(); + hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B"); + + if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded; + if(type == kFlagZ) fh1FFRatioZFolded[i] = hRatioUnfolded; + if(type == kFlagXi) fh1FFRatioXiFolded[i] = hRatioUnfolded; + + + // ratio backfolded to original histo + TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded); + hRatioBackFolded->Reset(); + hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B"); + + if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded; + if(type == kFlagZ) fh1FFRatioZBackFolded[i] = hRatioBackFolded; + if(type == kFlagXi) fh1FFRatioXiBackFolded[i] = hRatioBackFolded; + + delete hnHist; + delete hnFlatEfficiency; + + } +} + +//_____________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors) +{ + + Int_t type = kFlagPt; + UnfoldHistos(nIter, useCorrelatedErrors, type); +} + +//_____________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors) +{ + + Int_t type = kFlagZ; + UnfoldHistos(nIter, useCorrelatedErrors, type); +} + +//_____________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors) +{ + + Int_t type = kFlagXi; + UnfoldHistos(nIter, useCorrelatedErrors, type); +} + +//______________________________________________________________________________________________ +TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse) +{ + // apply (multiply) response matrix to hist + + TH1F* hOut = new TH1F(*hist); + hOut->Reset(); + + const Int_t axisM = 0; + const Int_t axisT = 1; + + Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins(); + Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins(); + + // response matrix normalization + // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign + // so for unfolding proper we leave it to CORRFW ... + + Double_t normFacResponse[nBinsT]; + + for(Int_t iT=1; iT<=nBinsT; iT++){ + + Double_t sumResp = 0; + + for(Int_t iM=1; iM<=nBinsM; iM++){ + + Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM); + Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT); + + Double_t binCoord[] = {coordM,coordT}; + + Long64_t binIndex = hnResponse->GetBin(binCoord); + + Double_t resp = hnResponse->GetBinContent(binIndex); + + sumResp += resp; + } + + normFacResponse[iT] = 0; + if(sumResp) normFacResponse[iT] = 1/sumResp; + } + + + + for(Int_t iM=1; iM<=nBinsM; iM++){ + + Double_t contM = 0; + + for(Int_t iT=1; iT<=nBinsT; iT++){ + + Double_t contT = hist->GetBinContent(iT); + + Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM); + Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT); + + Double_t binCoord[] = {coordM,coordT}; + + Long64_t binIndex = hnResponse->GetBin(binCoord); + + Double_t resp = hnResponse->GetBinContent(binIndex); + + contM += resp*normFacResponse[iT]*contT; + } + + hOut->SetBinContent(iM,contM); + } + + return hOut; +} + +//_______________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist) +{ + // read reconstruction efficiency from file + // argument strlist optional - read from directory strdir if not specified + + // temporary histos to hold histos from file + TH1F* hEffPt[fNJetPtSlices]; + TH1F* hEffZ[fNJetPtSlices]; + TH1F* hEffXi[fNJetPtSlices]; + + for(Int_t i=0; i0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + for(Int_t i=0; i (fJetPtSlices->At(i)); + Int_t jetPtUpLim = static_cast (fJetPtSlices->At(i+1)); + + TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim)); + + + if(list){ + hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); + hEffZ[i] = (TH1F*) list->FindObject(strNameEffZ); + hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); + } + else{ + hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); + hEffZ[i] = (TH1F*) gDirectory->Get(strNameEffZ); + hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); + } + + if(!hEffPt[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data()); + } + + if(!hEffZ[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data()); + } + + if(!hEffXi[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data()); + } + + + if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hEffZ[i] = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray()); + + if(hEffPt[i]) hEffPt[i]->SetDirectory(0); + if(hEffZ[i]) hEffZ[i]->SetDirectory(0); + if(hEffXi[i]) hEffXi[i]->SetDirectory(0); + + } // jet slices loop + + f.Close(); + + for(Int_t i=0; i0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + for(Int_t i=0; i (fJetPtSlices->At(i)); + Int_t jetPtUpLim = static_cast (fJetPtSlices->At(i+1)); + + TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim)); + TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim)); + TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim)); + + + if(list){ + hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); + hEffZBgr[i] = (TH1F*) list->FindObject(strNameEffZBgr); + hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); + } + else{ + hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); + hEffZBgr[i] = (TH1F*) gDirectory->Get(strNameEffZBgr); + hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); + } + + if(!hEffPtBgr[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data()); + } + + if(!hEffZBgr[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data()); + } + + if(!hEffXiBgr[i]){ + Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data()); + } + + + if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hEffZBgr[i] = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray()); + + if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); + if(hEffZBgr[i]) hEffZBgr[i]->SetDirectory(0); + if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); + + } // jet slices loop + + f.Close(); + + for(Int_t i=0; iGetTrackPt(i); + TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i); + TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i); + + TString histNamePt = histPt->GetName(); + TString histNameZ = histZ->GetName(); + TString histNameXi = histXi->GetName(); + + + TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt); + hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,""); + + TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ); + hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,""); + + TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi); + hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,""); + + fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr); + } +} + +//___________________________________________________ +void AliFragmentationFunctionCorrections::EffCorrBgr() +{ + // apply efficiency correction to bgr distributions + + AddCorrectionLevelBgr("eff"); + + Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1); + + for(Int_t i=0; iGetTrackPt(i); + TH1F* histZ = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i); + TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i); + + TString histNamePt = histPt->GetName(); + TString histNameZ = histZ->GetName(); + TString histNameXi = histXi->GetName(); + + + TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt); + hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,""); + + TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ); + hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,""); + + TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi); + hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,""); + + fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr); + } +} + +//______________________________________________________________________ +void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel) +{ + // re-evaluate jet energy after FF corrections from dN/dpt distribution + // apply correction (shift) to dN/dxi distribution: xi = ln (pt/E) -> xi' = ln (pt/E') = ln (pt/E x E/E') = xi + ln E/E' + // argument corrlevel: which level of correction to be corrected/shifted to + + if(corrLevel>=fNCorrectionLevels){ + Printf(" calc xi shift: corrLevel exceeded - do nothing"); + return; + } + + Double_t* jetPtUncorr = new Double_t[fNJetPtSlices]; + Double_t* jetPtCorr = new Double_t[fNJetPtSlices]; + Double_t* deltaXi = new Double_t[fNJetPtSlices]; + + for(Int_t i=0; iGetTrackPt(i); + TH1F* histPt = fCorrFF[corrLevel]->GetTrackPt(i); + + Double_t ptUncorr = 0; + Double_t ptCorr = 0; + + for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){ + + Double_t cont = histPtRaw->GetBinContent(bin); + Double_t width = histPtRaw->GetBinWidth(bin); + Double_t meanPt = histPtRaw->GetBinCenter(bin); + + ptUncorr += meanPt*cont*width; + } + + for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){ + + Double_t cont = histPt->GetBinContent(bin); + Double_t width = histPt->GetBinWidth(bin); + Double_t meanPt = histPt->GetBinCenter(bin); + + ptCorr += meanPt*cont*width; + } + + jetPtUncorr[i] = ptUncorr; + jetPtCorr[i] = ptCorr; + } + + // calc dXi from dN/dpt distribution : + // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim); + + Double_t ptUncorr = jetPtUncorr[i]; + Double_t ptCorr = jetPtCorr[i]; + + Double_t dXi = TMath::Log(ptCorr/ptUncorr); + + Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f " + ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi); + + deltaXi[i] = dXi; + } + + // book & fill new dN/dxi histos + + for(Int_t i=0; iGetXi(i); + + Double_t dXi = deltaXi[i]; + + Int_t nBins = histXi->GetNbinsX(); + const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray(); + Float_t binsVecNew[nBins+1]; + + TString strName = histXi->GetName(); + strName.Append("_shift"); + TString strTit = histXi->GetTitle(); + + TH1F* hXiCorr; + + // create shifted histo ... + + if(binsVec){ // binsVec only neq NULL if histo was rebinned before + + for(Int_t bin=0; binGetXaxis()->GetXmin(); + Double_t xMax = histXi->GetXaxis()->GetXmax(); + + xMin += dXi; + xMax += dXi; + + hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax); + } + + // ... and fill + + for(Int_t bin=1; binGetBinContent(bin); + Double_t err = histXi->GetBinError(bin); + + hXiCorr->SetBinContent(bin,cont); + hXiCorr->SetBinError(bin,err); + } + + new(fh1FFXiShift[i]) TH1F(*hXiCorr); + delete hXiCorr; + } + + delete[] jetPtUncorr; + delete[] jetPtCorr; + delete[] deltaXi; +} + +//_____________________________________________________ +void AliFragmentationFunctionCorrections::SubtractBgr() +{ + // subtract bgr distribution from FF + // the current corr level is used for both + + AddCorrectionLevel("bgrSub"); + + for(Int_t i=0; iGetTrackPt(i); + TH1F* histZ = fCorrFF[fNCorrectionLevels-2]->GetZ(i); + TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i); + + TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i); + TH1F* histZBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i); + TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i); + + TString histNamePt = histPt->GetName(); + TString histNameZ = histZ->GetName(); + TString histNameXi = histXi->GetName(); + + + TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt); + hFFTrackPtBgrSub->Add(histPtBgr,-1); + + TH1F* hFFZBgrSub = (TH1F*) histZ->Clone(histNameZ); + hFFZBgrSub->Add(histZBgr,-1); + + TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi); + hFFXiBgrSub->Add(histXiBgr,-1); + + fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub); + } +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile, TString strOutDir,TString strPostfix) +{ + // read task ouput from MC and write single track eff - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix); +} + +//___________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix) +{ + // read task output from MC and write single track eff as function of pt, eta, phi + // argument strlist optional - read from directory strdir if not specified + // write eff to file stroutfile - by default only 'update' file (don't overwrite) + + + TH1D* hdNdptTracksMCPrimGen; + TH2D* hdNdetadphiTracksMCPrimGen; + + TH1D* hdNdptTracksMCPrimRec; + TH2D* hdNdetadphiTracksMCPrimRec; + + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()){ + if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data()); + gDirectory->cd(strdir); + } + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + + TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen"; + if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); + + TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec"; + if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data()); + + TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen"; + if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data()); + + TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec"; + if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data()); + + + if(list){ + hdNdptTracksMCPrimGen = (TH1D*) list->FindObject(hnamePtRecEffGen); + hdNdetadphiTracksMCPrimGen = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen); + + hdNdptTracksMCPrimRec = (TH1D*) list->FindObject(hnamePtRecEffRec); + hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec); + } + else{ + hdNdptTracksMCPrimGen = (TH1D*) gDirectory->Get(hnamePtRecEffGen); + hdNdetadphiTracksMCPrimGen = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen); + + hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get(hnamePtRecEffRec); + hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec); + } + + hdNdptTracksMCPrimGen->SetDirectory(0); + hdNdetadphiTracksMCPrimGen->SetDirectory(0); + hdNdptTracksMCPrimRec->SetDirectory(0); + hdNdetadphiTracksMCPrimRec->SetDirectory(0); + + f.Close(); + + // projections: dN/deta, dN/dphi + + TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen"); + TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen"); + + TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec"); + TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec"); + + // rebin + + TString strNamePtGen = "hTrackPtGenPrim"; + TString strNamePtRec = "hTrackPtRecPrim"; + + if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray()); + if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray()); + + // efficiency: divide + + TString hNameTrackEffPt = "hSingleTrackEffPt"; + if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data()); + + TString hNameTrackEffEta = "hSingleTrackEffEta"; + if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data()); + + TString hNameTrackEffPhi = "hSingleTrackEffPhi"; + if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data()); + + + TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt); + hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors + + TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta); + hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors + + TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi); + hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors + + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + + if(strOutDir && strOutDir.Length()){ + + TDirectory* dir; + if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); + else{ + dir = out.mkdir(strOutDir); + dir->cd(); + } + } + + hSingleTrackEffPt->Write(); + hSingleTrackEffEta->Write(); + hSingleTrackEffPhi->Write(); + + out.Close(); +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile, TString strOutDir) +{ + // read task ouput from MC and write single track eff - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir); +} + +//___________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile, TString strOutDir) +{ + // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi + // argument strlist optional - read from directory strdir if not specified + // write corr factor to file stroutfile - by default only 'update' file (don't overwrite) + + TH1D* hdNdptTracksMCPrimRec; + TH2D* hdNdetadphiTracksMCPrimRec; + + TH1D* hdNdptTracksMCSecRec; + TH2D* hdNdetadphiTracksMCSecRec; + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + + if(list){ + hdNdptTracksMCPrimRec = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen"); + hdNdetadphiTracksMCPrimRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen"); + + hdNdptTracksMCSecRec = (TH1D*) list->FindObject("fh1TrackQAPtSecRec"); + hdNdetadphiTracksMCSecRec = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec"); + } + else{ + hdNdptTracksMCPrimRec = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen"); + hdNdetadphiTracksMCPrimRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen"); + + hdNdptTracksMCSecRec = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec"); + hdNdetadphiTracksMCSecRec = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec"); + } + + hdNdptTracksMCPrimRec->SetDirectory(0); + hdNdetadphiTracksMCPrimRec->SetDirectory(0); + hdNdptTracksMCSecRec->SetDirectory(0); + hdNdetadphiTracksMCSecRec->SetDirectory(0); + + f.Close(); + + // projections: dN/deta, dN/dphi + + TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec"); + TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec"); + + TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec"); + TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec"); + + + // rebin + + TString strNamePtPrim = "hTrackPtPrim"; + TString strNamePtSec = "hTrackPtSec"; + + if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray()); + if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray()); + + + // secondary correction factor: divide prim/(prim+sec) + + TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt"); + TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt"); + hSumPrimSecPt->Add(hdNdptTracksMCPrimRec); + hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors + + TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta"); + TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta"); + hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec); + hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors + + TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi"); + TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi"); + hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec); + hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors + + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + + if(strOutDir && strOutDir.Length()){ + + TDirectory* dir; + if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); + else{ + dir = out.mkdir(strOutDir); + dir->cd(); + } + } + + hdNdptTracksMCSecRec->Write(); + hdNdetaTracksMCSecRec->Write(); + hdNdphiTracksMCSecRec->Write(); + + hSingleTrackSecCorrPt->Write(); + hSingleTrackSecCorrEta->Write(); + hSingleTrackSecCorrPhi->Write(); + + out.Close(); +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile, TString strOutDir) +{ + // read task ouput from MC and write single track eff - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir); +} + + +//_____________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile, TString strOutDir) +{ + // read 2d THnSparse response matrices in pt from file + // project TH2 + // write to strOutfile + + THnSparse* hnResponseSinglePt; + TH2F* h2ResponseSinglePt; + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt"); + else hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt"); + + + if(!hnResponseSinglePt){ + Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__); + return; + } + + f.Close(); + + + // project 2d histo + h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim + h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); + + + // write + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + + if(strOutDir && strOutDir.Length()){ + + TDirectory* dir; + if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); + else{ + dir = out.mkdir(strOutDir); + dir->cd(); + } + } + + hnResponseSinglePt->Write(); + h2ResponseSinglePt->Write(); + + out.Close(); +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile) +{ + // read task ouput from MC and write single track eff - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile); +} + +//___________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile) +{ + // read task output from MC and write track eff in jet pt slices as function of pt, z, xi + // argument strlist optional - read from directory strdir if not specified + // write eff to file strOutfile - by default only 'update' file (don't overwrite) + + TH1F* hEffPt[fNJetPtSlices]; + TH1F* hEffXi[fNJetPtSlices]; + TH1F* hEffZ[fNJetPtSlices]; + + TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices]; + TH1F* hdNdxiMCPrimGen[fNJetPtSlices]; + TH1F* hdNdzMCPrimGen[fNJetPtSlices]; + + TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices]; + TH1F* hdNdxiMCPrimRec[fNJetPtSlices]; + TH1F* hdNdzMCPrimRec[fNJetPtSlices]; + + + TH1F* fh1FFJetPtRecEffGen; + + TH2F* fh2FFTrackPtRecEffGen; + TH2F* fh2FFZRecEffGen; + TH2F* fh2FFXiRecEffGen; + + TH2F* fh2FFTrackPtRecEffRec; + TH2F* fh2FFZRecEffRec; + TH2F* fh2FFXiRecEffRec; + + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + if(list){ + fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen"); + + fh2FFTrackPtRecEffGen = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen"); + fh2FFZRecEffGen = (TH2F*) list->FindObject("fh2FFZRecEffGen"); + fh2FFXiRecEffGen = (TH2F*) list->FindObject("fh2FFXiRecEffGen"); + + fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec"); + fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec"); + fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec"); + } + else{ + fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen"); + + fh2FFTrackPtRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen"); + fh2FFZRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen"); + fh2FFXiRecEffGen = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen"); + + fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec"); + fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec"); + fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec"); + } + + fh1FFJetPtRecEffGen->SetDirectory(0); + + fh2FFTrackPtRecEffGen->SetDirectory(0); + fh2FFZRecEffGen->SetDirectory(0); + fh2FFXiRecEffGen->SetDirectory(0); + + fh2FFTrackPtRecEffRec->SetDirectory(0); + fh2FFZRecEffRec->SetDirectory(0); + fh2FFXiRecEffRec->SetDirectory(0); + + f.Close(); + + + // projections: FF for generated and reconstructed primaries + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1; + + if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + // project + // appendix 'unbinned' to avoid histos with same name after rebinning + + hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + hdNdzMCPrimGen[i] = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o"); + hdNdxiMCPrimGen[i] = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o"); + + hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o"); + hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o"); + + // rebin + + if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hdNdzMCPrimGen[i] = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray()); + + if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray()); + + hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,""); + hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,""); + hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,""); + + hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,""); + hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,""); + hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,""); + + // normalize + + Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp); + + NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); + NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); + NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); + + NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); + NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); + NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); + + // divide rec/gen : efficiency + + TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + + hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt); + hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors + + hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi); + hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors + + hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ); + hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors + } + + // write + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + +// if(strdir && strdir.Length()){ +// TDirectory* dir = out.mkdir(strdir); +// dir->cd(); +// } + + for(Int_t i=0; iWrite(); + hdNdxiMCPrimGen[i]->Write(); + hdNdzMCPrimGen[i]->Write(); + + hdNdptTracksMCPrimRec[i]->Write(); + hdNdxiMCPrimRec[i]->Write(); + hdNdzMCPrimRec[i]->Write(); + + hEffPt[i]->Write(); + hEffXi[i]->Write(); + hEffZ[i]->Write(); + } + + out.Close(); + +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile) +{ + // read task ouput from MC and write secondary correction - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile); +} + +//___________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile) +{ + // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi + // argument strlist optional - read from directory strdir if not specified + // write eff to file strOutfile - by default only 'update' file (don't overwrite) + + TH1F* hSecCorrPt[fNJetPtSlices]; + TH1F* hSecCorrXi[fNJetPtSlices]; + TH1F* hSecCorrZ[fNJetPtSlices]; + + TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices]; + TH1F* hdNdxiMCPrimRec[fNJetPtSlices]; + TH1F* hdNdzMCPrimRec[fNJetPtSlices]; + + TH1F* hdNdptTracksMCSecRec[fNJetPtSlices]; + TH1F* hdNdxiMCSecRec[fNJetPtSlices]; + TH1F* hdNdzMCSecRec[fNJetPtSlices]; + + TH1F* fh1FFJetPtRecEffGen; + + TH2F* fh2FFTrackPtRecEffRec; + TH2F* fh2FFZRecEffRec; + TH2F* fh2FFXiRecEffRec; + + TH2F* fh2FFTrackPtSecRec; + TH2F* fh2FFZSecRec; + TH2F* fh2FFXiSecRec; + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + if(list){ + fh1FFJetPtRecEffGen = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen"); + + fh2FFTrackPtRecEffRec = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec"); + fh2FFZRecEffRec = (TH2F*) list->FindObject("fh2FFZRecEffRec"); + fh2FFXiRecEffRec = (TH2F*) list->FindObject("fh2FFXiRecEffRec"); + + fh2FFTrackPtSecRec = (TH2F*) list->FindObject("fh2FFTrackPtSecRec"); + fh2FFZSecRec = (TH2F*) list->FindObject("fh2FFZSecRec"); + fh2FFXiSecRec = (TH2F*) list->FindObject("fh2FFXiSecRec"); + } + else{ + fh1FFJetPtRecEffGen = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen"); + + fh2FFTrackPtRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec"); + fh2FFZRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec"); + fh2FFXiRecEffRec = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec"); + + fh2FFTrackPtSecRec = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec"); + fh2FFZSecRec = (TH2F*) gDirectory->FindObject("fh2FFZSecRec"); + fh2FFXiSecRec = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec"); + } + + fh1FFJetPtRecEffGen->SetDirectory(0); + + fh2FFTrackPtRecEffRec->SetDirectory(0); + fh2FFZRecEffRec->SetDirectory(0); + fh2FFXiRecEffRec->SetDirectory(0); + + fh2FFTrackPtSecRec->SetDirectory(0); + fh2FFZSecRec->SetDirectory(0); + fh2FFXiSecRec->SetDirectory(0); + + f.Close(); + + + // projections: FF for generated and reconstructed primaries + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1; + + if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + // project + // appendix 'unbinned' to avoid histos with same name after rebinning + + hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + hdNdzMCPrimRec[i] = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o"); + hdNdxiMCPrimRec[i] = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o"); + + hdNdptTracksMCSecRec[i] = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range + hdNdzMCSecRec[i] = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o"); + hdNdxiMCSecRec[i] = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o"); + + // rebin + + if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hdNdzMCPrimRec[i] = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray()); + + if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hdNdzMCSecRec[i] = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i] = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray()); + + hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,""); + hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,""); + hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,""); + + hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,""); + hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,""); + hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,""); + + // normalize + + Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp); + + NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); + NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); + NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); + + NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin); + NormalizeTH1(hdNdzMCSecRec[i],nJetsBin); + NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin); + + // divide rec/gen : efficiency + + TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + + hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt); + TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt"); + hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]); + hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors + + hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi); + TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi"); + hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]); + hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors + + hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ); + TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ"); + hSumPrimSecZ->Add(hdNdzMCPrimRec[i]); + hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors + } + + // write + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + + for(Int_t i=0; iWrite(); + // hdNdxiMCSecRec[i]->Write(); + // hdNdzMCSecRec[i]->Write(); + + hSecCorrPt[i]->Write(); + hSecCorrXi[i]->Write(); + hSecCorrZ[i]->Write(); + } + + out.Close(); +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile, + Bool_t updateOutfile) +{ + // read task ouput from MC and write single track eff - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile); +} + +//_____________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist, + TString strOutfile, Bool_t updateOutfile) +{ + // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file + // project THnSparse + TH2 in jet pt slices + // write to strOutfile + + THnSparse* hn3ResponseJetPt; + THnSparse* hn3ResponseJetZ; + THnSparse* hn3ResponseJetXi; + + // 2D response matrices + + THnSparse* hnResponsePt[fNJetPtSlices]; + THnSparse* hnResponseZ[fNJetPtSlices]; + THnSparse* hnResponseXi[fNJetPtSlices]; + + TH2F* h2ResponsePt[fNJetPtSlices]; + TH2F* h2ResponseZ[fNJetPtSlices]; + TH2F* h2ResponseXi[fNJetPtSlices]; + + // 1D projections on gen pt / rec pt axes + + TH1F* h1FFPtRec[fNJetPtSlices]; + TH1F* h1FFZRec[fNJetPtSlices]; + TH1F* h1FFXiRec[fNJetPtSlices]; + + TH1F* h1FFPtGen[fNJetPtSlices]; + TH1F* h1FFZGen[fNJetPtSlices]; + TH1F* h1FFXiGen[fNJetPtSlices]; + + TH1F* h1RatioPt[fNJetPtSlices]; + TH1F* h1RatioZ[fNJetPtSlices]; + TH1F* h1RatioXi[fNJetPtSlices]; + + + + TFile f(strInfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + if(list){ + hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt"); + hn3ResponseJetZ = (THnSparse*) list->FindObject("fhnResponseJetZ"); + hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi"); + } + else{ + hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt"); + hn3ResponseJetZ = (THnSparse*) gDirectory->Get("fhnResponseJetZ"); + hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi"); + } + + + if(!hn3ResponseJetPt){ + Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__); + return; + } + + if(!hn3ResponseJetZ){ + Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__); + return; + } + + if(!hn3ResponseJetXi){ + Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__); + return; + } + + f.Close(); + + // axes + + Int_t axisJetPtTHn3 = -1; + Int_t axisGenPtTHn3 = -1; + Int_t axisRecPtTHn3 = -1; + + for(Int_t i=0; iGetNdimensions(); i++){ + + TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle(); + + if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i; + if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i; + if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i; + } + + if(axisJetPtTHn3 == -1){ + Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__); + return; + } + + if(axisGenPtTHn3 == -1){ + Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__); + return; + } + + if(axisRecPtTHn3 == -1){ + Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__); + return; + } + + for(Int_t i=0; iAt(i); + Float_t jetPtUpLim = fJetPtSlices->At(i+1); + + Int_t binLo = static_cast(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim)); + Int_t binUp = static_cast(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1; + + if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){ + Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); + return; + } + + TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast (jetPtLoLim),static_cast (jetPtUpLim))); + + TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + + TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + + TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast(jetPtLoLim),static_cast(jetPtUpLim))); + + + // 2D projections in jet pt range + + hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); + hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); + hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); + + Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3}; + + hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj); + hnResponseZ[i] = hn3ResponseJetZ->Projection(2,axesProj); + hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj); + + hnResponsePt[i]->SetNameTitle(strNameRespPt,""); + hnResponseZ[i]->SetNameTitle(strNameRespZ,""); + hnResponseXi[i]->SetNameTitle(strNameRespXi,""); + + h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim + h2ResponseZ[i] = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim + h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim + + h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,""); + h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,""); + h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,""); + + + // 1D projections + + Int_t axisGenPtTHn2 = -1; + Int_t axisRecPtTHn2 = -1; + + for(Int_t d=0; dGetNdimensions(); d++){ + + TString title = hnResponsePt[i]->GetAxis(d)->GetTitle(); + + if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d; + if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d; + } + + + if(axisGenPtTHn2 == -1){ + Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__); + return; + } + + if(axisRecPtTHn2 == -1){ + Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__); + return; + } + + + h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim + h1FFZRec[i] = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim + h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim + + h1FFPtRec[i]->SetNameTitle(strNameRecPt,""); + h1FFZRec[i]->SetNameTitle(strNameRecZ,""); + h1FFXiRec[i]->SetNameTitle(strNameRecXi,""); + + h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim + h1FFZGen[i] = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim + h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim + + h1FFPtGen[i]->SetNameTitle(strNameGenPt,""); + h1FFZGen[i]->SetNameTitle(strNameGenZ,""); + h1FFXiGen[i]->SetNameTitle(strNameGenXi,""); + + // normalize 1D projections + + if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) h1FFZRec[i] = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray()); + + if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) h1FFZGen[i] = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray()); + + NormalizeTH1(h1FFPtRec[i],fNJets->At(i)); + NormalizeTH1(h1FFZRec[i],fNJets->At(i)); + NormalizeTH1(h1FFXiRec[i],fNJets->At(i)); + + NormalizeTH1(h1FFPtGen[i],fNJets->At(i)); + NormalizeTH1(h1FFZGen[i],fNJets->At(i)); + NormalizeTH1(h1FFXiGen[i],fNJets->At(i)); + + // ratios 1D projections + + h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt); + h1RatioPt[i]->Reset(); + h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B"); + + h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ); + h1RatioZ[i]->Reset(); + h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B"); + + h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi); + h1RatioXi[i]->Reset(); + h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B"); + } + + + // write + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data()); + + //if(strdir && strdir.Length()){ + // TDirectory* dir = out.mkdir(strdir); + // dir->cd(); + //} + + for(Int_t i=0; iWrite(); + hnResponseXi[i]->Write(); + hnResponseZ[i]->Write(); + + h2ResponsePt[i]->Write(); + h2ResponseXi[i]->Write(); + h2ResponseZ[i]->Write(); + + h1FFPtRec[i]->Write(); + h1FFZRec[i]->Write(); + h1FFXiRec[i]->Write(); + + h1FFPtGen[i]->Write(); + h1FFZGen[i]->Write(); + h1FFXiGen[i]->Write(); + + h1RatioPt[i]->Write(); + h1RatioZ[i]->Write(); + h1RatioXi[i]->Write(); + + } + + out.Close(); +} + +//______________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist) +{ + // read response matrices from file + // argument strlist optional - read from directory strdir if not specified + // note: THnSparse are not rebinned + + THnSparse* hRespPt[fNJetPtSlices]; + THnSparse* hRespZ[fNJetPtSlices]; + THnSparse* hRespXi[fNJetPtSlices]; + + for(Int_t i=0; i0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + for(Int_t i=0; i (fJetPtSlices->At(i)); + Int_t jetPtUpLim = static_cast (fJetPtSlices->At(i+1)); + + TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim)); + + if(list){ + hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt); + hRespZ[i] = (THnSparse*) list->FindObject(strNameRespZ); + hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi); + } + else{ + hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt); + hRespZ[i] = (THnSparse*) gDirectory->Get(strNameRespZ); + hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi); + } + + if(!hRespPt[i]){ + Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data()); + } + + if(!hRespZ[i]){ + Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data()); + } + + if(!hRespXi[i]){ + Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data()); + } + + // if(0){ // can't rebin THnSparse ... + // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray()); + // if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray()); + + // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray()); + // if(fNHistoBinsZ[i]) hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray()); + + // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray()); + // if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray()); + // } + + + } // jet slices loop + + f.Close(); // THnSparse pointers still valid even if file closed + +// for(Int_t i=0; i0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + // temporary histos to store pointers from file + TH1F* hist[fNJetPtSlices]; + + for(Int_t i=0; i (fJetPtSlices->At(i)); + Int_t jetPtUpLim = static_cast (fJetPtSlices->At(i+1)); + + TString strName; + + if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); + if(type == kFlagZ) strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); + if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); + + hist[i] = (TH1F*) gDirectory->Get(strName); + + if(!hist[i]){ + Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data()); + } + + + //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray()); + + if(hist[i]) hist[i]->SetDirectory(0); + + } // jet slices loop + + f.Close(); + + + for(Int_t i=0; iGetTrackPt(i); // levels -1: latest corr level +// TH1F* histZRec = fCorrFF[fNCorrectionLevels-1]->GetZ(i); // levels -1: latest corr level +// TH1F* histXiRec = fCorrFF[fNCorrectionLevels-1]->GetXi(i); // levels -1: latest corr level + +// TH1F* histPtGen = fh1FFTrackPtGenPrim[i]; +// TH1F* histZGen = fh1FFZGenPrim[i]; +// TH1F* histXiGen = fh1FFXiGenPrim[i]; + +// TString histNamePt = histPtRec->GetName(); +// TString histNameZ = histZRec->GetName(); +// TString histNameXi = histXiRec->GetName(); + +// histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen"); +// histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen"); +// histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen"); + +// // ratio +// TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt); +// hRatioRecGenPt->Reset(); +// hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B"); + +// TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ); +// hRatioRecGenZ->Reset(); +// hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B"); + +// TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi); +// hRatioRecGenXi->Reset(); +// hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B"); + +// new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt); +// new(fh1FFRatioRecGenZ[i]) TH1F(*hRatioRecGenZ); +// new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi); +// } +// } + +// //___________________________________________________________ +// void AliFragmentationFunctionCorrections::RatioRecPrimaries() +// { +// // create ratio reconstructed tracks over reconstructed primaries +// // use raw FF (corrLevel 0) + +// Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__); + +// for(Int_t i=0; iGetTrackPt(i); // levels -1: latest corr level +// TH1F* histZRec = fCorrFF[corrLevel]->GetZ(i); // levels -1: latest corr level +// TH1F* histXiRec = fCorrFF[corrLevel]->GetXi(i); // levels -1: latest corr level + +// TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i]; +// TH1F* histZRecPrim = fh1FFZRecPrim[i]; +// TH1F* histXiRecPrim = fh1FFXiRecPrim[i]; + +// TString histNamePt = histPtRec->GetName(); +// TString histNameZ = histZRec->GetName(); +// TString histNameXi = histXiRec->GetName(); + +// histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim"); +// histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim"); +// histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim"); + +// // ratio +// TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt); +// hRatioRecPrimPt->Reset(); +// hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B"); + +// TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ); +// hRatioRecPrimZ->Reset(); +// hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B"); + +// TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi); +// hRatioRecPrimXi->Reset(); +// hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B"); + + +// new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt); +// new(fh1FFRatioRecPrimZ[i]) TH1F(*hRatioRecPrimZ); +// new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi); +// } +// } + +// __________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile) +{ + + // project response matrices on both axes: + // FF for rec primaries, in terms of generated and reconstructed momentum + // write FF and ratios to outFile + + Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data()); + + TH1F* hFFPtRec[fNJetPtSlices]; + TH1F* hFFZRec[fNJetPtSlices]; + TH1F* hFFXiRec[fNJetPtSlices]; + + TH1F* hFFPtGen[fNJetPtSlices]; + TH1F* hFFZGen[fNJetPtSlices]; + TH1F* hFFXiGen[fNJetPtSlices]; + + TH1F* hRatioPt[fNJetPtSlices]; + TH1F* hRatioZ[fNJetPtSlices]; + TH1F* hRatioXi[fNJetPtSlices]; + + + Int_t axisGenPt = 1; + Int_t axisRecPt = 0; + + for(Int_t i=0; i (fJetPtSlices->At(i)); + Int_t jetPtUpLim = static_cast (fJetPtSlices->At(i+1)); + + TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + + TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim)); + + TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim)); + TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim)); + + + hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim + hFFZRec[i] = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim + hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim + + hFFPtRec[i]->SetNameTitle(strNameRecPt,""); + hFFZRec[i]->SetNameTitle(strNameRecZ,""); + hFFXiRec[i]->SetNameTitle(strNameRecXi,""); + + + hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim + hFFZGen[i] = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim + hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim + + hFFPtGen[i]->SetNameTitle(strNameGenPt,""); + hFFZGen[i]->SetNameTitle(strNameGenZ,""); + hFFXiGen[i]->SetNameTitle(strNameGenXi,""); + + + if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hFFZRec[i] = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hFFXiRec[i] = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray()); + + if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); + if(fNHistoBinsZ[i]) hFFZGen[i] = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray()); + if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray()); + + NormalizeTH1(hFFPtGen[i],fNJets->At(i)); + NormalizeTH1(hFFZGen[i],fNJets->At(i)); + NormalizeTH1(hFFXiGen[i],fNJets->At(i)); + + NormalizeTH1(hFFPtRec[i],fNJets->At(i)); + NormalizeTH1(hFFZRec[i],fNJets->At(i)); + NormalizeTH1(hFFXiRec[i],fNJets->At(i)); + + + hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt); + hRatioPt[i]->Reset(); + hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B"); + + hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ); + hRatioZ[i]->Reset(); + hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B"); + + hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi); + hRatioXi[i]->Reset(); + hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B"); + } + + + + // write + + TFile out(strOutfile,"RECREATE"); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + for(Int_t i=0; iWrite(); + hFFZRec[i]->Write(); + hFFXiRec[i]->Write(); + + hFFPtGen[i]->Write(); + hFFZGen[i]->Write(); + hFFXiGen[i]->Write(); + + hRatioPt[i]->Write(); + hRatioZ[i]->Write(); + hRatioXi[i]->Write(); + } + + out.Close(); +} + +// ____________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir ) +{ + // project response matrix on both axes: + // pt spec for rec primaries, in terms of generated and reconstructed momentum + // write spec and ratios to outFile + + Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data()); + + TH1F* hSpecPtRec; + TH1F* hSpecPtGen; + TH1F* hRatioPt; + + Int_t axisGenPt = 1; + Int_t axisRecPt = 0; + + TString strNameRecPt = "h1SpecTrackPtRecPrim_recPt"; + TString strNameGenPt = "h1SpecTrackPtRecPrim_genPt"; + TString strNameRatioPt = "h1RatioTrackPtRecPrim"; + + hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim + hSpecPtRec->SetNameTitle(strNameRecPt,""); + + hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim + hSpecPtGen->SetNameTitle(strNameGenPt,""); + + hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt); + hRatioPt->Reset(); + hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B"); + + TString outfileOption = "RECREATE"; + if(updateOutfile) outfileOption = "UPDATE"; + + TFile out(strOutfile,outfileOption); + + if(!out.IsOpen()){ + Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data()); + return; + } + + + if(strOutDir && strOutDir.Length()){ + + TDirectory* dir; + if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); + else{ + dir = out.mkdir(strOutDir); + dir->cd(); + } + } + + hSpecPtRec->Write(); + hSpecPtGen->Write(); + hRatioPt->Write(); + + out.Close(); +} + + +//__________________________________________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type) +{ + // rebin histo, rescale bins according to new width + // only correct for input histos with equal bin size + + // args: jetPtSlice, type, use current corr level + + // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...) + // array size of binsLimits: nBinsLimits + // array size of binsWidth: nBinsLimits-1 + // binsLimits have to be in increasing order + // if binning undefined for any slice, original binning will be kept + + if(!fNJetPtSlices){ + Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__); + return; + } + + if(jetPtSlice>=fNJetPtSlices){ + Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice); + return; + } + + + Double_t binLimitMin = binsLimits[0]; + Double_t binLimitMax = binsLimits[nBinsLimits-1]; + + Double_t binLimit = binLimitMin; // start value + + Int_t sizeUpperLim = 1000; //static_cast(binLimitMax/binsWidth[0])+1; - only if first bin has smallest width, but not case for dN/dxi ... + TArrayD binsArray(sizeUpperLim); + Int_t nBins = 0; + binsArray.SetAt(binLimitMin,nBins++); + + while(binLimit= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor + } + + Double_t currentBinWidth = binsWidth[currentSlice]; + binLimit += currentBinWidth; + + binsArray.SetAt(binLimit,nBins++); + } + + + TH1F* hist = 0; + if(type == kFlagPt) hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice); + if(type == kFlagZ) hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice); + if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice); + if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0); + + + Double_t binWidthNoRebin = hist->GetBinWidth(1); + + Double_t* bins = binsArray.GetArray(); + + hist = (TH1F*) hist->Rebin(nBins-1,"",bins); + + for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){ + + Double_t binWidthRebin = hist->GetBinWidth(bin); + Double_t scaleF = binWidthNoRebin / binWidthRebin; + + Double_t binCont = hist->GetBinContent(bin); + Double_t binErr = hist->GetBinError(bin); + + binCont *= scaleF; + binErr *= scaleF; + + hist->SetBinContent(bin,binCont); + hist->SetBinError(bin,binErr); + } + + + + TH1F* temp = new TH1F(*hist); + + if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0); + if(type == kFlagZ) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0); + if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp); + if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0); + + + delete temp; +} +//__________________________________________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist, TString strOutfile) +{ + + if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + THnSparse* hn6ResponseJetPt; + + if(list){ + hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation"); + } + else{ + hn6ResponseJetPt = (THnSparse*) list->FindObject("fhnCorrelation"); + } + + Int_t axis6RecJetPt = 0; + Int_t axis6GenJetPt = 3; + + hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)"); + hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)"); + + Int_t nBinsRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins(); + Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1); + Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt); + + Int_t nBinsGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins(); + Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1); + Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt); + + Int_t nBinsTrackPt = 200; + Int_t loLimTrackPt = 0; + Int_t upLimTrackPt = 200; + + + Int_t nBinsResponse[4] = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt}; + Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt}; + Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt}; + + const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"}; + + THnSparseD* hn4ResponseTrackPtJetPt = new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse); + + for(Int_t i=0; i<4; i++){ + hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]); + } + + + // fill + + +} + +//_____________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname) +{ + + ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency); + +} + +//_____________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname) +{ + + ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse); + +} + +//_____________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname) +{ + + ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries); + +} + +//______________________________________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type) +{ + // read single track correction (pt) from file + // type: efficiency / response / secondaries correction + + if(!((type == kFlagEfficiency) || (type == kFlagResponse) || (type == kFlagSecondaries))){ + Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__); + return; + } + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0 && type==kFlagEfficiency) Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + if(fDebug>0 && type==kFlagResponse) Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + if(strdir && strdir.Length()) gDirectory->cd(strdir); + + TList* list = 0; + + if(strlist && strlist.Length()){ + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + } + + TH1F* h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later + THnSparse* hnCorrHist = 0; + + if(type == kFlagEfficiency || type == kFlagSecondaries){ + + if(list) h1CorrHist = (TH1F*) list->FindObject(strname); + else h1CorrHist = (TH1F*) gDirectory->Get(strname); + + if(!h1CorrHist){ + Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data()); + return; + } + + } + else if(type == kFlagResponse){ + + if(list) hnCorrHist = (THnSparse*) list->FindObject(strname); + else hnCorrHist = (THnSparse*) gDirectory->Get(strname); + + if(!hnCorrHist){ + Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data()); + return; + } + + } + + if(h1CorrHist) h1CorrHist->SetDirectory(0); + //if(hnCorrHist) hnCorrHist->SetDirectory(0); + + f.Close(); + + if(type == kFlagEfficiency) fh1EffSinglePt = h1CorrHist; + else if(type == kFlagResponse) fhnResponseSinglePt = hnCorrHist; + else if(type == kFlagSecondaries) fh1SecCorrSinglePt = h1CorrHist; + +} + +//________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID) +{ + // read track pt spec from task ouput - standard dir/list + + TString strdir = "PWG4_FragmentationFunction_" + strID; + TString strlist = "fracfunc_" + strID; + + ReadRawPtSpec(strInfile,strdir,strlist); +} + +//_______________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist) +{ + // get raw pt spectra from file + + + // book histos + fNCorrectionLevelsSinglePt = 0; + fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels]; + AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum + + // get raw pt spec from input file, normalize + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + gDirectory->cd(strdir); + + TList* list = 0; + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + + TString hnameTrackPt("fh1TrackQAPtRecCuts"); + TString hnameEvtSel("fh1EvtSelection"); + + TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt); + TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel); + + if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; } + if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; } + + + //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0)); + + // evts after physics selection + Float_t nEvents = fh1EvtSel->GetEntries() + - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1)) // evts rejected by trigger selection ( = PhysicsSelection + - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(2)); // evts with wrong centrality class + + + fh1TrackPt->SetDirectory(0); + + f.Close(); + + + NormalizeTH1(fh1TrackPt,nEvents); + + // raw FF = corr level 0 + fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt); +} + + +//_______________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist) +{ + // get raw pt spectra from file + // for output from Martas QA task + + + // book histos + fNCorrectionLevelsSinglePt = 0; + fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels]; + AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum + + // get raw pt spec from input file, normalize + + TFile f(strfile,"READ"); + + if(!f.IsOpen()){ + Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data()); + return; + } + + if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data()); + + gDirectory->cd(strdir); + + TList* list = 0; + + if(!(list = (TList*) gDirectory->Get(strlist))){ + Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data()); + return; + } + + TString hnameTrackPt("fPtSel"); + TString hnameEvtSel("fNEventAll"); + + TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt); + TH1F* fh1EvtSel = (TH1F*) list->FindObject(hnameEvtSel); + + if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; } + if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; } + + + // evts after physics selection + Float_t nEvents = fh1EvtSel->GetEntries(); + + fh1TrackPt->SetDirectory(0); + + f.Close(); + + + NormalizeTH1(fh1TrackPt,nEvents); + + // raw FF = corr level 0 + fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt); +} + +// ________________________________________________________ +void AliFragmentationFunctionCorrections::EffCorrSinglePt() +{ + // apply efficiency correction to inclusive track pt spec + + AddCorrectionLevelSinglePt("eff"); + + TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); + + if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__); + + TString histNamePt = histPt->GetName(); + TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt); + hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,""); + + fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr); +} + +//___________________________________________________________________________________________________________________________ +void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors) +{ + // unfolde inclusive dN/dpt spectra + + AddCorrectionLevelSinglePt("unfold"); + + TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded + THnSparse* hnResponse = fhnResponseSinglePt; + + TString histNameTHn = hist->GetName(); + if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn"); + if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt"); + + + TString histNameBackFolded = hist->GetName(); + histNameBackFolded.Append("_backfold"); + + TString histNameRatioFolded = hist->GetName(); + if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio"); + if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt"); + histNameRatioFolded.Append("_unfold"); + + TString histNameRatioBackFolded = hist->GetName(); + if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio"); + if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt"); + histNameRatioBackFolded.Append("_backfold"); + + THnSparse* hnHist = TH1toSparse(hist,histNameTHn,hist->GetTitle()); + THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff + + THnSparse* hnUnfolded + = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors); + + TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); + hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle()); + + fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded); + + // backfolding: apply response matrix to unfolded spectrum + TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); + hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle()); + + fh1SingleTrackPtBackFolded = hBackFolded; + + + // ratio unfolded to original histo + TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded); + hRatioUnfolded->Reset(); + hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B"); + + fh1RatioSingleTrackPtFolded = hRatioUnfolded; + + + // ratio backfolded to original histo + TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded); + hRatioBackFolded->Reset(); + hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B"); + + fh1RatioSingleTrackPtBackFolded = hRatioBackFolded; + + delete hnHist; + delete hnFlatEfficiency; + +} + +// ________________________________________________________ +void AliFragmentationFunctionCorrections::SecCorrSinglePt() +{ + // apply efficiency correction to inclusive track pt spec + + AddCorrectionLevelSinglePt("secCorr"); + + TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); + + if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX()) + Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__); + + TString histNamePt = histPt->GetName(); + TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt); + + hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,""); + + fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr); +} diff --git a/PWG4/JetTasks/AliFragmentationFunctionCorrections.h b/PWG4/JetTasks/AliFragmentationFunctionCorrections.h new file mode 100644 index 00000000000..d66eea80f35 --- /dev/null +++ b/PWG4/JetTasks/AliFragmentationFunctionCorrections.h @@ -0,0 +1,242 @@ +// ***************************************************************************** +// * Task for corrections to output from AliAnalysisTaskFragmentationFunctions * +// **************************************************************************** + +#ifndef ALIFRAGMENTATIONFUNCTIONCORRECTIONS_H +#define ALIFRAGMENTATIONFUNCTIONCORRECTIONS_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +#include "TObject.h" + +class ThnSparse; + +class AliFragmentationFunctionCorrections : public TObject { + + public: + + //---------------------------------------- + class AliFragFuncCorrHistos : public TObject + { + + public: + + AliFragFuncCorrHistos(); + AliFragFuncCorrHistos(const AliFragFuncCorrHistos& copy); + AliFragFuncCorrHistos& operator=(const AliFragFuncCorrHistos &o); + virtual ~AliFragFuncCorrHistos(); + AliFragFuncCorrHistos(const char* label,Int_t arraySize); + void AddCorrHistos(Int_t slice,TH1F* histPt=0,TH1F* histZ=0,TH1F* histXi=0); + void ReplaceCorrHistos(Int_t slice,TH1F* histPt=0,TH1F* histZ=0,TH1F* histXi=0); + + TH1F* GetTrackPt(const Int_t slice); + TH1F* GetZ(const Int_t slice); + TH1F* GetXi(const Int_t slice); + + TString GetLabel() { return fCorrLabel; } + + private: + + Int_t fArraySize; + + TH1F** fh1CorrFFTrackPt; //! corrected FF histos + TH1F** fh1CorrFFZ; //! corrected FF histos + TH1F** fh1CorrFFXi; //! corrected FF histos + + TString fCorrLabel; //! correction label + + ClassDef(AliFragFuncCorrHistos, 1); + }; + + AliFragmentationFunctionCorrections(); + AliFragmentationFunctionCorrections(const AliFragmentationFunctionCorrections ©); + AliFragmentationFunctionCorrections& operator=(const AliFragmentationFunctionCorrections &o); + virtual ~AliFragmentationFunctionCorrections(); + + virtual void SetDebugLevel(Int_t debug){ fDebug = debug; } + + void DeleteHistoArray(TH1F** hist) const; + void DeleteTHnSparseArray(THnSparse** hist) const; + TH1F** BookHistoArray(); + THnSparse** BookTHnSparseArray(); + void AddCorrectionLevel(const char* label = ""); + void AddCorrectionLevelBgr(const char* label = ""); + void AddCorrectionLevelSinglePt(const char* label = ""); + + void SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices); + + void SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins,Int_t type); + void SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth,Int_t type); + void SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins); + void SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth); + + // set histo bins for inclusive pt spectra + + void NormalizeTH1(TH1* hist, const Float_t nJets); + void NormalizeFF(); + void NormalizeBgr(); + void ReadRawFF(TString strfile, TString strID, TString strFFID = "RecCuts"); + void ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID); + void ReadRawBgr(TString strfile, TString strID, TString strBgrID = "Perp", TString strFFID = "RecCuts"); + void ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID); + void ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID); + void ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID); + + void WriteOutput(TString strfile, TString strdir = "", Bool_t updateOutfile = kTRUE); + + THnSparse* TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst = kFALSE); + + THnSparse* Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter, + const Bool_t useCorrelatedErrors = kTRUE, const THnSparse* hnPrior = 0x0); + + void UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type); + + void UnfoldPt(const Int_t nIter=5, const Bool_t useCorrelatedErrors=kTRUE); + void UnfoldZ(const Int_t nIter=5, const Bool_t useCorrelatedErrors=kTRUE); + void UnfoldXi(const Int_t nIter=5, const Bool_t useCorrelatedErrors=kTRUE); + + TH1F* ApplyResponse(const TH1F* hist, THnSparse* hnResponse); + + void ReadEfficiency(TString strfile, TString strdir = "", TString strlist = ""); + void ReadBgrEfficiency(TString strfile, TString strdir = "", TString strlist = ""); + + void EffCorr(); + void EffCorrBgr(); + + void XiShift(const Int_t corrLevel); + + void SubtractBgr(); + + void WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE, TString strOutDir = "", TString strPostfix = ""); + void WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, TString strOutfile, Bool_t updateOutfile = kTRUE, TString strOutDir = "", + TString strPostfix = ""); + + void WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE, TString strOutDir = ""); + void WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, TString strOutfile, Bool_t updateOutfile = kTRUE, TString strOutDir = ""); + + void WriteSingleResponse(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE, TString strOutDir = ""); + void WriteSingleResponse(TString strInfile, TString strdir, TString strlist, TString strOutfile, Bool_t updateOutfile = kTRUE, TString strOutDir = ""); + + void WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE); + void WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, TString strOutfile, Bool_t updateOutfile = kTRUE); + + void WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE); + void WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, TString strOutfile, Bool_t updateOutfile = kTRUE); + + void WriteJetResponse(TString strInfile, TString strID, TString strOutfile,Bool_t updateOutfile = kTRUE); + void WriteJetResponse(TString strInfile, TString strdir, TString strlist,TString strOutfile, Bool_t updateOutfile); + + void ReadResponse(TString strfile, TString strdir="", TString strlist=""); + void ReadPriors(TString strfile,const Int_t type); + + void ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir = ""); + void ProjectJetResponseMatrices(TString strOutfile); + + void RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type); + + void WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist, TString strOutfile); + + void ReadSingleTrackEfficiency(TString strfile, TString strdir="", TString strlist="", TString strname="hSingleTrackEffPt"); + void ReadSingleTrackResponse(TString strfile, TString strdir="", TString strlist="", TString strname="fhnResponseSinglePt"); + void ReadSingleTrackSecCorr(TString strfile, TString strdir="", TString strlist="", TString strname="hSingleTrackSecCorrPt"); + void ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type); + + void ReadRawPtSpec(TString strInfile, TString strID); + void ReadRawPtSpec(TString strfile, TString strdir, TString strlist); + void ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist); // spectra from Martas QA task + void EffCorrSinglePt(); + void UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors); + void SecCorrSinglePt(); + + + enum {kFlagPt=0,kFlagZ,kFlagXi,kFlagSinglePt}; + enum {kFlagEfficiency=0,kFlagResponse,kFlagSecondaries}; + + + private: + + static const Int_t fgMaxNCorrectionLevels = 10; //! max number of corrections + + Int_t fDebug; //! Debug level + Int_t fNJetPtSlices; //! n slices in jet pt + TArrayF* fJetPtSlices; //! array to hold slices in jet pt + + TArrayF* fNJets; //! nJets per jet pt slice - non-int e.g. for xsec/nTrials scaled spectra + TArrayF* fNJetsBgr; //! nJets bgr per jet pt slice - non-int e.g. for xsec/nTrials scaled spectra + + Int_t fNHistoBinsSinglePt; //! nBins inclusive pt spec histos - leave undefinded to use original binning + TArrayD* fHistoBinsSinglePt; //! inclusive pt spec histo bins + + Int_t* fNHistoBinsPt; //! nBins FF histos in each jet pt slice - leave undefinded for any slice to use original binning + Int_t* fNHistoBinsZ; //! nBins FF histos in each jet pt slice - leave undefinded for any slice to use original binning + Int_t* fNHistoBinsXi; //! nBins FF histos in each jet pt slice - leave undefinded for any slice to use original binning + + TArrayD** fHistoBinsPt; //! FF histo bins + TArrayD** fHistoBinsZ; //! FF histo bins + TArrayD** fHistoBinsXi; //! FF histo bins + + Int_t fNCorrectionLevels; //! corrections applied: efficiency, secondaries, resolution unfolding, bgr subtraction + AliFragFuncCorrHistos** fCorrFF; //! array of fragmentation functions, dimensions: jet pt bins, correction steps + + Int_t fNCorrectionLevelsBgr; //! corrections applied: efficiency, secondaries, resolution unfolding, bgr subtraction + AliFragFuncCorrHistos** fCorrBgr; //! array of bgr fragmentation functions, dimensions: jet pt bins, correction steps + + Int_t fNCorrectionLevelsSinglePt; //! corrections applied: efficiency, secondaries, resolution unfolding, bgr subtraction + AliFragFuncCorrHistos** fCorrSinglePt; //! array to keep single track pt spectra, 1D in jet pt bins dimension + + + + // xi shift + TH1F** fh1FFXiShift; //! FF: track xi, corrected for shift in jet energy + + // eff correction + TH1F* fh1EffSinglePt; //! efficiency all tracks + + TH1F** fh1EffPt; //! reconstruction efficiency track pt + TH1F** fh1EffZ; //! reconstruction efficiency z + TH1F** fh1EffXi; //! reconstruction efficiency xi + + TH1F** fh1EffBgrPt; //! reconstruction efficiency bgr track pt + TH1F** fh1EffBgrZ; //! reconstruction efficiency bgr z + TH1F** fh1EffBgrXi; //! reconstruction efficiency bgr xi + + + // unfolding + + TH1F** fh1FFTrackPtBackFolded; //! FF: track pt backfolded (unfolded + smeared with response matrix) + TH1F** fh1FFZBackFolded; //! FF: track z, backfolded (unfolded + smeared with response matrix) + TH1F** fh1FFXiBackFolded; //! FF: track xi,backfolded (unfolded + smeared with response matrix) + + TH1F** fh1FFRatioTrackPtFolded; //! ratio FF: track pt unfolded / original + TH1F** fh1FFRatioZFolded; //! ratio FF: track z unfolded / original + TH1F** fh1FFRatioXiFolded; //! ratio FF: track xi unfolded / original + + TH1F** fh1FFRatioTrackPtBackFolded; //! ratio FF: track pt backfolded / original + TH1F** fh1FFRatioZBackFolded; //! ratio FF: track z backfolded / original + TH1F** fh1FFRatioXiBackFolded; //! ratio FF: track xi backfolded / original + + TH1F* fh1SingleTrackPtBackFolded; //! inclusive track pt backfolded (unfolded + smeared with response matrix) + TH1F* fh1RatioSingleTrackPtFolded; //! ratio inclusive track pt unfolded / original + TH1F* fh1RatioSingleTrackPtBackFolded; //! ratio inblusive track pt backfolded / original + + THnSparse* fhnResponseSinglePt; //! response matrix pt gen vs rec all tracks + THnSparse** fhnResponsePt; //! response matrix pt gen vs rec + THnSparse** fhnResponseZ; //! response matrix z gen vs rec + THnSparse** fhnResponseXi; //! response matrix xi gen vs rec + + TH1F** fh1FFTrackPtPrior; //! FF: track pt prior + TH1F** fh1FFZPrior; //! FF: track z prior + TH1F** fh1FFXiPrior; //! FF: track xi prior + + // secondaries + TH1F* fh1SecCorrSinglePt; //! secondaries correction all tracks + + + + ClassDef(AliFragmentationFunctionCorrections, 1); +}; + +#endif diff --git a/PWG4/PWG4JetTasksLinkDef.h b/PWG4/PWG4JetTasksLinkDef.h index b5db74e61af..97279fbeb5b 100644 --- a/PWG4/PWG4JetTasksLinkDef.h +++ b/PWG4/PWG4JetTasksLinkDef.h @@ -27,6 +27,8 @@ #pragma link C++ class AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos+; #pragma link C++ class AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos+; #pragma link C++ class AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos+; +#pragma link C++ class AliFragmentationFunctionCorrections+; +#pragma link C++ class AliFragmentationFunctionCorrections::AliFragFuncCorrHistos+; #pragma link C++ class AliAnalysisTaskJetChem+; #pragma link C++ class AliAnalysisTaskJetChem::AliFragFuncHistosInvMass+; #pragma link C++ class AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass+;