X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FAliJetResponseMaker.cxx;h=b8731048df988fe2dc7cac64baeec3d584b56f62;hb=a976bee95f39765a266649587aef594f540804b8;hp=1d1dd1823791196f6a2d9c74539039166dcf6e5c;hpb=3c9775d958e21b0bd985e052e70c6c1d519e365f;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx index 1d1dd182379..b8731048df9 100644 --- a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx +++ b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx @@ -37,18 +37,12 @@ AliJetResponseMaker::AliJetResponseMaker() : fDeltaEtaDeltaPhiAxis(0), fNEFAxis(0), fZAxis(0), - fDoJet2Histogram(0), fIsJet1Rho(kFALSE), fIsJet2Rho(kFALSE), - fHistLeadingJets1PtArea(0), - fHistLeadingJets1CorrPtArea(0), - fHistLeadingJets2PtArea(0), - fHistLeadingJets2CorrPtArea(0), - fHistLeadingJets2PtAreaAcceptance(0), - fHistLeadingJets2CorrPtAreaAcceptance(0), + fHistRejectionReason1(0), + fHistRejectionReason2(0), fHistJets1(0), fHistJets2(0), - fHistJets2Acceptance(0), fHistMatching(0), fHistJets1PhiEta(0), fHistJets1PtArea(0), @@ -59,9 +53,6 @@ AliJetResponseMaker::AliJetResponseMaker() : fHistJets2PhiEta(0), fHistJets2PtArea(0), fHistJets2CorrPtArea(0), - fHistJets2PhiEtaAcceptance(0), - fHistJets2PtAreaAcceptance(0), - fHistJets2CorrPtAreaAcceptance(0), fHistJets2NEFvsPt(0), fHistJets2CEFvsCEFPt(0), fHistJets2ZvsPt(0), @@ -127,18 +118,12 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) : fDeltaEtaDeltaPhiAxis(0), fNEFAxis(0), fZAxis(0), - fDoJet2Histogram(0), fIsJet1Rho(kFALSE), fIsJet2Rho(kFALSE), - fHistLeadingJets1PtArea(0), - fHistLeadingJets1CorrPtArea(0), - fHistLeadingJets2PtArea(0), - fHistLeadingJets2CorrPtArea(0), - fHistLeadingJets2PtAreaAcceptance(0), - fHistLeadingJets2CorrPtAreaAcceptance(0), + fHistRejectionReason1(0), + fHistRejectionReason2(0), fHistJets1(0), fHistJets2(0), - fHistJets2Acceptance(0), fHistMatching(0), fHistJets1PhiEta(0), fHistJets1PtArea(0), @@ -149,9 +134,6 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) : fHistJets2PhiEta(0), fHistJets2PtArea(0), fHistJets2CorrPtArea(0), - fHistJets2PhiEtaAcceptance(0), - fHistJets2PtAreaAcceptance(0), - fHistJets2CorrPtAreaAcceptance(0), fHistJets2NEFvsPt(0), fHistJets2CEFvsCEFPt(0), fHistJets2ZvsPt(0), @@ -276,27 +258,6 @@ void AliJetResponseMaker::AllocateTH2() fOutput->Add(fHistJets2CorrPtArea); } - fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2); - fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta"); - fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistJets2PhiEtaAcceptance); - - fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", - fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt); - fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area"); - fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); - fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistJets2PtAreaAcceptance); - - if (fIsJet2Rho) { - fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", - fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt); - fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area"); - fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)"); - fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts"); - fOutput->Add(fHistJets2CorrPtAreaAcceptance); - } - fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJets2ZvsPt->GetXaxis()->SetTitle("Z"); fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)"); @@ -635,10 +596,10 @@ void AliJetResponseMaker::AllocateTHnSparse() { // Allocate THnSparse histograms. - TString title[20]= {""}; - Int_t nbins[20] = {0}; - Double_t min[20] = {0.}; - Double_t max[20] = {0.}; + TString title[25]= {""}; + Int_t nbins[25] = {0}; + Double_t min[25] = {0.}; + Double_t max[25] = {0.}; Int_t dim = 0; title[dim] = "#phi"; @@ -665,16 +626,26 @@ void AliJetResponseMaker::AllocateTHnSparse() max[dim] = 1.5; dim++; - title[dim] = "NEF"; - nbins[dim] = fNbins/4; - min[dim] = 0; - max[dim] = 1.2; - dim++; + if (fNEFAxis) { + title[dim] = "NEF"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + } - title[dim] = "Z"; - nbins[dim] = fNbins/4; + if (fZAxis) { + title[dim] = "Z"; + nbins[dim] = fNbins/4; + min[dim] = 0; + max[dim] = 1.2; + dim++; + } + + title[dim] = "p_{T,particle}^{leading} (GeV/c)"; + nbins[dim] = 120; min[dim] = 0; - max[dim] = 1.2; + max[dim] = 120; dim++; Int_t dim1 = dim, dim2 = dim; @@ -708,17 +679,10 @@ void AliJetResponseMaker::AllocateTHnSparse() dim2++; } - if (fDoJet2Histogram) { - fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max); - for (Int_t i = 0; i < dim2; i++) - fHistJets2->GetAxis(i)->SetTitle(title[i]); - fOutput->Add(fHistJets2); - } - - fHistJets2Acceptance = new THnSparseD("fHistJets2Acceptance","fHistJets2Acceptance",dim2,nbins,min,max); + fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max); for (Int_t i = 0; i < dim2; i++) - fHistJets2Acceptance->GetAxis(i)->SetTitle(title[i]); - fOutput->Add(fHistJets2Acceptance); + fHistJets2->GetAxis(i)->SetTitle(title[i]); + fOutput->Add(fHistJets2); // Matching @@ -766,6 +730,18 @@ void AliJetResponseMaker::AllocateTHnSparse() max[dim] = 1.2; dim++; + title[dim] = "p_{T,particle,1}^{leading} (GeV/c)"; + nbins[dim] = 120; + min[dim] = 0; + max[dim] = 120; + dim++; + + title[dim] = "p_{T,particle,2}^{leading} (GeV/c)"; + nbins[dim] = 120; + min[dim] = 0; + max[dim] = 120; + dim++; + if (fDeltaPtAxis) { title[dim] = "#deltaA_{jet}"; nbins[dim] = fNbins/2; @@ -882,6 +858,20 @@ void AliJetResponseMaker::UserCreateOutputObjects() if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE; else fIsJet2Rho = kTRUE; + fHistRejectionReason1 = new TH2F("fHistRejectionReason1", "fHistRejectionReason1", 32, 0, 32, 100, 0, 250); + fHistRejectionReason1->GetXaxis()->SetTitle("Rejection reason"); + fHistRejectionReason1->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); + fHistRejectionReason1->GetZaxis()->SetTitle("counts"); + SetRejectionReasonLabels(fHistRejectionReason1->GetXaxis()); + fOutput->Add(fHistRejectionReason1); + + fHistRejectionReason2 = new TH2F("fHistRejectionReason2", "fHistRejectionReason2", 32, 0, 32, 100, 0, 250); + fHistRejectionReason2->GetXaxis()->SetTitle("Rejection reason"); + fHistRejectionReason2->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)"); + fHistRejectionReason2->GetZaxis()->SetTitle("counts"); + SetRejectionReasonLabels(fHistRejectionReason2->GetXaxis()); + fOutput->Add(fHistRejectionReason2); + if (fHistoType==0) AllocateTH2(); else @@ -891,7 +881,7 @@ void AliJetResponseMaker::UserCreateOutputObjects() } //________________________________________________________________________ -void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set) +void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set) { if (fHistoType==1) { THnSparse *histo = 0; @@ -899,8 +889,6 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, histo = fHistJets1; else if (Set==2) histo = fHistJets2; - else if (Set==3) - histo = fHistJets2Acceptance; if (!histo) return; @@ -920,11 +908,13 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, else if (title=="NEF") contents[i] = NEF; else if (title=="Z") - contents[i] = Z; + contents[i] = LeadingPt/Pt; else if (title=="p_{T}^{corr}") contents[i] = CorrPt; else if (title=="p_{T}^{MC}") contents[i] = MCPt; + else if (title=="p_{T,particle}^{leading} (GeV/c)") + contents[i] = LeadingPt; else AliWarning(Form("Unable to fill dimension %s!",title.Data())); } @@ -935,7 +925,7 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, if (Set == 1) { fHistJets1PtArea->Fill(A, Pt); fHistJets1PhiEta->Fill(Eta, Phi); - fHistJets1ZvsPt->Fill(Z, Pt); + fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt); fHistJets1NEFvsPt->Fill(NEF, Pt); fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); if (fIsJet1Rho) @@ -944,17 +934,11 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, else if (Set == 2) { fHistJets2PtArea->Fill(A, Pt); fHistJets2PhiEta->Fill(Eta, Phi); - if (fIsJet2Rho) - fHistJets2CorrPtArea->Fill(A, CorrPt); - } - else if (Set == 3) { - fHistJets2PtAreaAcceptance->Fill(A, Pt); - fHistJets2PhiEtaAcceptance->Fill(Eta, Phi); - fHistJets2ZvsPt->Fill(Z, Pt); + fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt); fHistJets2NEFvsPt->Fill(NEF, Pt); fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt); if (fIsJet2Rho) - fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt); + fHistJets2CorrPtArea->Fill(A, CorrPt); } } } @@ -962,7 +946,7 @@ void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, //________________________________________________________________________ void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, - Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t Z1, Double_t Z2) + Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2) { if (fHistoType==1) { Double_t contents[20]={0}; @@ -1006,9 +990,13 @@ void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_ else if (title=="NEF_{2}") contents[i] = NEF2; else if (title=="Z_{1}") - contents[i] = Z1; + contents[i] = LeadingPt1/Pt1; else if (title=="Z_{2}") - contents[i] = Z2; + contents[i] = LeadingPt2/Pt2; + else if (title=="p_{T,particle,1}^{leading} (GeV/c)") + contents[i] = LeadingPt1; + else if (title=="p_{T,particle,2}^{leading} (GeV/c)") + contents[i] = LeadingPt2; else AliWarning(Form("Unable to fill dimension %s!",title.Data())); } @@ -1222,7 +1210,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet } Int_t MClabel = TMath::Abs(track->GetLabel()); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely @@ -1248,7 +1236,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely @@ -1269,7 +1257,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet clus->GetMomentum(part, const_cast(fVertex)); Int_t MClabel = TMath::Abs(clus->GetLabel()); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel != 0) continue; // this is not a MC particle; remove it completely @@ -1291,7 +1279,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet continue; } Int_t MClabel = TMath::Abs(track->GetLabel()); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel <= 0) continue; Int_t index = -1; @@ -1332,7 +1320,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell); Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId)); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel <= 0) continue; Int_t index1 = -1; @@ -1369,7 +1357,7 @@ void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet clus->GetMomentum(part, const_cast(fVertex)); Int_t MClabel = TMath::Abs(clus->GetLabel()); - if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift; + MClabel -= fMCLabelShift; if (MClabel <= 0) continue; Int_t index = -1; @@ -1644,15 +1632,19 @@ Bool_t AliJetResponseMaker::FillHistograms() jets2->ResetCurrentID(); while ((jet2 = jets2->GetNextJet())) { + AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID())); - if (fDoJet2Histogram) - FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), - jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 2); + if (jet2->Pt() < jets2->GetJetPtCut()) continue; + + Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2); + Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area(); if (jets2->AcceptJet(jet2)) - FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), - jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 3); + FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2, + corrpt2, jet2->MCPt(), 2); + else + fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(), jet2->Pt()); jet1 = jet2->MatchedJet(); @@ -1676,21 +1668,28 @@ Bool_t AliJetResponseMaker::FillHistograms() ce2 = jet2->ClosestJetDistance(); } + Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); - Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area(); FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), - jet1->NEF(), jet2->NEF(), jet1->MaxPartPt()/jet1->Pt(), jet2->MaxPartPt()/jet2->Pt()); + jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2); } jets1->ResetCurrentID(); - while ((jet1 = jets1->GetNextAcceptJet())) { + while ((jet1 = jets1->GetNextJet())) { + if (!jets1->AcceptJet(jet1)) { + fHistRejectionReason1->Fill(jets1->GetRejectionReasonBitPosition(), jet1->Pt()); + continue; + } if (jet1->MCPt() < fMinJetMCPt) continue; AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID())); - FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), jet1->MaxPartPt()/jet1->Pt(), - jet1->Pt() - jets1->GetRhoVal() * jet1->Area(), jet1->MCPt(), 1); + Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1); + Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area(); + + FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1, + corrpt1, jet1->MCPt(), 1); } return kTRUE; }