diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
index afd95b8b64faa10fac8a067f0fc3a201259c6a9e..e5fb3add05a44a4ff47345b5ade8ec7e8dd452be 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
@@ -69,6 +69,85 @@
#include "L1Field.h"
+
+/*------------------------------------------------------- clang 11.0.0 version (crashed)
+=======
+>>>>>>> Update combinatorial BG procedure
+
+#include "CbmAnaDielectronTask.h"
+#include "CbmGlobalTrack.h"
+#include "CbmKF.h"
+
+<<<<<<< HEAD
+#include "CbmL1PFFitter.h"
+
+=======
+>>>>>>> Update combinatorial BG procedure
+#include "CbmLmvmCandidate.h"
+#include "CbmLmvmUtils.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMvdHit.h"
+#include "CbmRichHit.h"
+#include "CbmRichPoint.h"
+#include "CbmRichRing.h"
+<<<<<<< HEAD
+
+#include "CbmL1PFFitter.h"
+#include "CbmStsHit.h"
+
+=======
+#include "CbmL1PFFitter.h"
+#include "CbmStsHit.h"
+>>>>>>> Update combinatorial BG procedure
+#include "CbmStsKFTrackFitter.h"
+#include "CbmStsTrack.h"
+#include "CbmTofHit.h"
+#include "CbmTofPoint.h"
+#include "CbmTrackMatchNew.h"
+#include "CbmTrdHit.h"
+#include "CbmTrdTrack.h"
+#include "CbmVertex.h"
+#include "cbm/elid/CbmLitGlobalElectronId.h"
+
+#include "FairBaseParSet.h"
+#include "FairEventHeader.h"
+#include "FairGeoMedium.h"
+#include "FairGeoNode.h"
+#include "FairGeoTransform.h"
+#include "FairGeoVector.h"
+#include "FairGeoVolume.h"
+#include "FairMCPoint.h"
+#include "FairRootManager.h"
+#include "FairRunAna.h"
+#include "FairRuntimeDb.h"
+#include "FairTask.h"
+#include "FairTrackParam.h"
+
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
+#include "TF1.h"
+#include "TGraph.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TObjArray.h"
+#include "TObject.h"
+#include "TProfile.h"
+#include "TRandom3.h"
+#include "TStopwatch.h"
+#include "TString.h"
+#include "TSystem.h"
+#include "TVector3.h"
+#include <TFile.h>
+
+#include <sstream>
+#include <vector>
+
+#include "L1Field.h"
+-----------------------------------------------------*/
+
using namespace std;
ClassImp(CbmAnaDielectronTask);
@@ -166,6 +245,7 @@ CbmAnaDielectronTask::CbmAnaDielectronTask()
, fUseTrd(kTRUE)
, fUseTof(kTRUE)
, fCandidates()
+ , fCandidatesTotal()
, fSTCandidates()
, fTTCandidates()
, fRTCandidates()
@@ -189,6 +269,16 @@ CbmAnaDielectronTask::CbmAnaDielectronTask()
, fh_vertex_el_gamma_rz()
, fh_signal_minv()
, fh_bg_minv()
+ , fh_combPairsPM_minv_sameEvent()
+ , fh_combPairsPP_minv_sameEvent()
+ , fh_combPairsMM_minv_sameEvent()
+ , fh_combPairsPM_minv_mixedEvents()
+ , fh_combPairsPP_minv_mixedEvents()
+ , fh_combPairsMM_minv_mixedEvents()
+ , fh_nof_plutoElectrons()
+ , fh_nof_plutoPositrons()
+ , fh_nof_urqmdElectrons()
+ , fh_nof_urqmdPositrons()
, fh_pi0_minv()
, fh_eta_minv()
, fh_gamma_minv()
@@ -235,6 +325,7 @@ CbmAnaDielectronTask::CbmAnaDielectronTask()
, fh_source_pairs_epem()
, fh_source_pairs(NULL)
, fh_event_number(NULL)
+ , fh_event_number_mixed(NULL)
, fh_nof_bg_tracks(NULL)
, fh_nof_el_tracks(NULL)
, fh_source_tracks(NULL)
@@ -379,6 +470,8 @@ void CbmAnaDielectronTask::InitHists()
// Event number counter
fh_event_number = new TH1D("fh_event_number", "fh_event_number", 1, 0, 1.);
fHistoList.push_back(fh_event_number);
+ fh_event_number_mixed = new TH1D("fh_event_number_mixed", "fh_event_number_mixed", 1, 0, 1.);
+ fHistoList.push_back(fh_event_number_mixed);
CreateSourceTypesH1(fh_richann, "fh_richann", "ANN output", "Yield", 100, -1.1, 1.1);
CreateSourceTypesH1(fh_trdann, "fh_trdann", "ANN output", "Yield", 100, -1.1, 1.1);
@@ -438,6 +531,31 @@ void CbmAnaDielectronTask::InitHists()
CreateAnalysisStepsH1(fh_pi0_minv, "fh_pi0_minv", "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0., 4.);
CreateAnalysisStepsH1(fh_eta_minv, "fh_eta_minv", "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0., 4.);
CreateAnalysisStepsH1(fh_gamma_minv, "fh_gamma_minv", "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0., 4.);
+
+ // Histograms for combinatorial BG
+ CreateAnalysisStepsH1(fh_combPairsPM_minv_sameEvent, "fh_combPairsPM_minv_sameEvent", "M_{e+e-} [GeV/c^{2}]", "Yield",
+ 4000, 0., 4.);
+ CreateAnalysisStepsH1(fh_combPairsPP_minv_sameEvent, "fh_combPairsPP_minv_sameEvent", "M_{e+e+} [GeV/c^{2}]", "Yield",
+ 4000, 0., 4.);
+ CreateAnalysisStepsH1(fh_combPairsMM_minv_sameEvent, "fh_combPairsMM_minv_sameEvent", "M_{e-e-} [GeV/c^{2}]", "Yield",
+ 4000, 0., 4.);
+
+ CreateAnalysisStepsH1(fh_combPairsPM_minv_mixedEvents, "fh_combPairsPM_minv_mixedEvents", "M_{e+e-} [GeV/c^{2}]",
+ "Yield", 4000, 0., 4.);
+ CreateAnalysisStepsH1(fh_combPairsPP_minv_mixedEvents, "fh_combPairsPP_minv_mixedEvents", "M_{e+e+} [GeV/c^{2}]",
+ "Yield", 4000, 0., 4.);
+ CreateAnalysisStepsH1(fh_combPairsMM_minv_mixedEvents, "fh_combPairsMM_minv_mixedEvents", "M_{e-e-} [GeV/c^{2}]",
+ "Yield", 4000, 0., 4.);
+
+ CreateAnalysisStepsH1(fh_nof_plutoElectrons, "fh_nof_plutoElectrons", "P [Gev/c]", "number pluto electrons", 10000, 0,
+ 10.);
+ CreateAnalysisStepsH1(fh_nof_plutoPositrons, "fh_nof_plutoPositrons", "P [Gev/c]", "number pluto positrons", 10000, 0,
+ 10.);
+ CreateAnalysisStepsH1(fh_nof_urqmdElectrons, "fh_nof_urqmdElectrons", "P [Gev/c]", "number urqmd electrons", 10000, 0,
+ 10.);
+ CreateAnalysisStepsH1(fh_nof_urqmdPositrons, "fh_nof_urqmdPositrons", "P [Gev/c]", "number urqmd positrons", 10000, 0,
+ 10.);
+
// minv for true matched and mismatched tracks
CreateAnalysisStepsH1(fh_bg_truematch_minv, "fh_bg_truematch_minv", "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0., 4.);
CreateAnalysisStepsH1(fh_bg_truematch_el_minv, "fh_bg_truematch_el_minv", "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0.,
@@ -453,8 +571,9 @@ void CbmAnaDielectronTask::InitHists()
CreateAnalysisStepsH1(fh_source_bg_minv[i], ss.str(), "M_{ee} [GeV/c^{2}]", "Yield", 4000, 0., 4.);
}
//Invariant mass vs. Mc Pt
- CreateAnalysisStepsH2(fh_signal_minv_pt, "fh_signal_minv_pt", "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", "Yield", 100, 0.,
- 2., 20, 0., 2.);
+ CreateAnalysisStepsH2(
+ fh_signal_minv_pt, "fh_signal_minv_pt", "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", "Yield", 100, 0., 4., 20, 0.,
+ 2.); // MIND: CHANGED "100, 0., 2." to 100, 0., 4., since in DrawAll.cxx this histo is added with the others below and compiler complains about various axes
CreateAnalysisStepsH2(fh_pi0_minv_pt, "fh_pi0_minv_pt", "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", "Yield", 100, 0., 4.,
20, 0., 2.);
CreateAnalysisStepsH2(fh_eta_minv_pt, "fh_eta_minv_pt", "M_{ee} [GeV/c^{2}]", "P_{t} [GeV/c]", "Yield", 100, 0., 4.,
@@ -666,6 +785,8 @@ void CbmAnaDielectronTask::Exec(Option_t*)
{
fh_event_number->Fill(0.5);
+ fEventNumber = fh_event_number->GetEntries();
+
Bool_t useMbias = false; // false for 40% central agag collisions (b<7.7fm)
Bool_t isCentralCollision = false;
@@ -675,7 +796,7 @@ void CbmAnaDielectronTask::Exec(Option_t*)
if (impactPar <= 7.7) isCentralCollision = true;
}
- cout << "-I- CbmAnaDielectronTask, event number " << fh_event_number->GetEntries() << endl;
+ cout << "-I- CbmAnaDielectronTask, event number " << fEventNumber << endl;
cout << "fPionMisidLevel = " << fPionMisidLevel << endl;
fCuts.Print();
cout << "fWeight = " << fWeight << endl;
@@ -700,7 +821,7 @@ void CbmAnaDielectronTask::Exec(Option_t*)
DifferenceSignalAndBg();
SignalAndBgReco();
FillElPiMomHist();
- FillNofChargedParticlesPerEvent();
+ FillNofChargedParticles();
FillMomLikeHist();
}
} // Exec
@@ -769,7 +890,7 @@ void CbmAnaDielectronTask::FillRichRingNofHits()
void CbmAnaDielectronTask::MCPairs()
{
- Int_t nMcTracks = fMCTracks->GetEntriesFast();
+ Int_t nMcTracks = fMCTracks->GetEntries();
for (Int_t i = 0; i < nMcTracks; i++) {
CbmMCTrack* mctrack = (CbmMCTrack*) fMCTracks->At(i);
Int_t motherId = mctrack->GetMotherId();
@@ -793,19 +914,6 @@ void CbmAnaDielectronTask::MCPairs()
fh_mc_signal_mom_angle->Fill(sqrtPMc, angle);
}
}
- for (Int_t iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
- if (i == iMc2) continue;
- CbmMCTrack* mctrack2 = (CbmMCTrack*) fMCTracks->At(iMc2);
- Int_t motherIdMc2 = mctrack2->GetMotherId();
- if (motherId == motherIdMc2 && CbmLmvmUtils::IsMcSignalElectron(mctrack2)) {
- CbmLmvmKinematicParams pKin = CbmLmvmKinematicParams::KinematicParamsWithMcTracks(mctrack, mctrack2);
- Double_t angle = pKin.fAngle;
- Double_t pMc = mctrack->GetP();
- Double_t pMc2 = mctrack2->GetP();
- Double_t sqrtPMc = TMath::Sqrt(pMc * pMc2);
- fh_mc_signal_mom_angle->Fill(sqrtPMc, angle);
- }
- }
}
if (isMcGammaElectron) {
TVector3 v;
@@ -867,28 +975,98 @@ void CbmAnaDielectronTask::RichPmtXY()
}
}
-void CbmAnaDielectronTask::FillNofChargedParticlesPerEvent()
+void CbmAnaDielectronTask::FillNofChargedParticles()
{
- Int_t nofMcTracks = fMCTracks->GetEntriesFast();
+ Int_t nofMcTracks = fMCTracks->GetEntries();
+ Int_t nCand = fCandidates.size();
Int_t nofChargedUrqmdParticles = 0;
Int_t nofChargedUrqmdParticlesAcc = 0;
+
+ cout << "FillNofChargedParticles: nofMcTracks = " << nofMcTracks << endl;
+ cout << "FillNofChargedParticles: nCand = " << nCand << endl;
+
for (Int_t i = 0; i < nofMcTracks; i++) {
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(i);
- Bool_t isMcTrackCharged = false;
- //if (mcTrack->GetCharge() != 0) isMcTrackCharged = true; // FIXME: TODO: uncomment when bug is fixed; issue https://lxcbmredmine01.gsi.de/issues/1826;
+ CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(i);
+ Int_t motherId = mcTrack->GetMotherId();
+ Int_t pdg = TMath::Abs(mcTrack->GetPdgCode());
+ bool isMcElectron = (pdg == 11) ? true : false;
+ double mom = mcTrack->GetP();
+ double charge =
+ mcTrack
+ ->GetCharge(); // FIXME: TODO: uncomment when bug is fixed; issue https://lxcbmredmine01.gsi.de/issues/1826; gives charge = +/- 3!
+ Bool_t isMcTrackCharged = (charge == 0) ? false : true;
Bool_t isMcSignalElectron = CbmLmvmUtils::IsMcSignalElectron(mcTrack);
- Bool_t isMcTrackAccepted = false;
- if (mcTrack->GetNPoints(ECbmModuleId::kSts) >= 4) isMcTrackAccepted = true;
- Bool_t isPrimary = false;
- Int_t motherId = mcTrack->GetMotherId();
- if (motherId == -1) isPrimary = true;
+ Bool_t isMcTrAcc =
+ IsMcTrackAccepted(i); // before: = (mcTrack->GetNPoints(ECbmModuleId::kSts) >= 4) ? true : false;
+ Bool_t isPrimary = (motherId == -1) ? true : false;
+
+ if (charge != 0) cout << "charge = " << charge << endl;
+
if (!isMcSignalElectron && isMcTrackCharged && isPrimary) nofChargedUrqmdParticles++;
- if (!isMcSignalElectron && isMcTrackCharged && isMcTrackAccepted && isPrimary) nofChargedUrqmdParticlesAcc++;
+ if (!isMcSignalElectron && isMcTrackCharged && isMcTrAcc && isPrimary) nofChargedUrqmdParticlesAcc++;
+
+ if (isMcSignalElectron && charge < 0)
+ fh_nof_plutoElectrons[kMc]->Fill(mom); // 'isMcEl' was redundant here and in next lines with cond. 'isMcSignalEl'
+ if (isMcSignalElectron && charge > 0) fh_nof_plutoPositrons[kMc]->Fill(mom);
+ if (isMcSignalElectron && isMcTrAcc && charge < 0) fh_nof_plutoElectrons[kAcc]->Fill(mom);
+ if (isMcSignalElectron && isMcTrAcc && charge > 0) fh_nof_plutoPositrons[kAcc]->Fill(mom);
+ if (!isMcSignalElectron && isMcElectron && charge < 0) fh_nof_urqmdElectrons[kMc]->Fill(mom);
+ if (!isMcSignalElectron && isMcElectron && charge > 0) fh_nof_urqmdPositrons[kMc]->Fill(mom);
+ if (!isMcSignalElectron && isMcElectron && isMcTrAcc && charge < 0) fh_nof_urqmdElectrons[kAcc]->Fill(mom);
+ if (!isMcSignalElectron && isMcElectron && isMcTrAcc && charge > 0) fh_nof_urqmdPositrons[kAcc]->Fill(mom);
}
fh_nof_charged_particles->Fill(nofChargedUrqmdParticles);
fh_nof_charged_particles_acc->Fill(nofChargedUrqmdParticlesAcc);
-}
+ for (Int_t i = 0; i < nCand; i++) {
+ int charge = fCandidates[i].fCharge;
+ Bool_t isMcSignalElectron = fCandidates[i].fIsMcSignalElectron;
+
+ Bool_t isChiPrimary = fCandidatesTotal[i].fChi2Prim < fCuts.fChiPrimCut;
+ Bool_t isElectron = fCandidatesTotal[i].fIsElectron;
+ Bool_t isGammaCut = !fCandidatesTotal[i].fIsGamma;
+ //Bool_t isMvd1Cut = fCandidatesTotal[i].fIsMvd1CutElectron;
+ //Bool_t isMvd2Cut = fCandidatesTotal[i].fIsMvd2CutElectron;
+ Bool_t isStCut = fCandidatesTotal[i].fIsStCutElectron;
+ Bool_t isRtCut = fCandidatesTotal[i].fIsRtCutElectron;
+ Bool_t isTtCut = fCandidatesTotal[i].fIsTtCutElectron;
+ Bool_t isPtCut = fCandidatesTotal[i].fMomentum.Perp() > fCuts.fPtCut;
+
+ if (isMcSignalElectron && charge < 0) fh_nof_plutoElectrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && charge > 0) fh_nof_plutoPositrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && charge < 0)
+ fh_nof_plutoElectrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && charge > 0)
+ fh_nof_plutoPositrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && charge < 0)
+ fh_nof_plutoElectrons[kTtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && charge > 0)
+ fh_nof_plutoPositrons[kTtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge < 0)
+ fh_nof_plutoElectrons[kPtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge > 0)
+ fh_nof_plutoPositrons[kPtCut]->Fill(fCandidates[i].fMomentum.Mag());
+
+ if (!isMcSignalElectron && charge < 0) fh_nof_urqmdElectrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && charge > 0) fh_nof_urqmdPositrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && charge < 0)
+ fh_nof_urqmdElectrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && charge > 0)
+ fh_nof_urqmdPositrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && charge < 0)
+ fh_nof_urqmdElectrons[kTtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && charge > 0)
+ fh_nof_urqmdPositrons[kTtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge < 0)
+ fh_nof_urqmdElectrons[kPtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ if (!isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge > 0)
+ fh_nof_urqmdPositrons[kPtCut]->Fill(fCandidates[i].fMomentum.Mag());
+ }
+}
Bool_t CbmAnaDielectronTask::IsMcTrackAccepted(Int_t mcTrackInd)
{
@@ -901,7 +1079,7 @@ Bool_t CbmAnaDielectronTask::IsMcTrackAccepted(Int_t mcTrackInd)
void CbmAnaDielectronTask::SingleParticleAcceptance()
{
- Int_t nMcTracks = fMCTracks->GetEntriesFast();
+ Int_t nMcTracks = fMCTracks->GetEntries();
for (Int_t i = 0; i < nMcTracks; i++) {
CbmMCTrack* mctrack = (CbmMCTrack*) fMCTracks->At(i);
Int_t motherId = mctrack->GetMotherId();
@@ -938,7 +1116,7 @@ void CbmAnaDielectronTask::SingleParticleAcceptance()
void CbmAnaDielectronTask::PairMcAndAcceptance()
{
- Int_t nMcTracks = fMCTracks->GetEntriesFast();
+ Int_t nMcTracks = fMCTracks->GetEntries();
for (Int_t iP = 0; iP < nMcTracks; iP++) {
CbmMCTrack* mctrackP = (CbmMCTrack*) fMCTracks->At(iP);
// Int_t motherIdP = mctrackP->GetMotherId();
@@ -989,7 +1167,7 @@ void CbmAnaDielectronTask::PairMcAndAcceptance()
void CbmAnaDielectronTask::FillElPiMomHist()
{
- Int_t nMcTracks = fMCTracks->GetEntriesFast();
+ Int_t nMcTracks = fMCTracks->GetEntries();
for (Int_t i = 0; i < nMcTracks; i++) {
CbmMCTrack* mctrack = (CbmMCTrack*) fMCTracks->At(i);
// Int_t motherId = mctrack->GetMotherId();
@@ -1147,6 +1325,7 @@ void CbmAnaDielectronTask::FillCandidates()
cand.fIsRtCutElectron = false;
cand.fIsMvd1CutElectron = true;
cand.fIsMvd2CutElectron = true;
+ cand.fEventNumber = fEventNumber;
CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iGTrack);
if (NULL == gTrack) continue;
@@ -1175,6 +1354,7 @@ void CbmAnaDielectronTask::FillCandidates()
if (pdg == 211) {
IsElectron(iGTrack, cand.fMomentum.Mag(), &cand);
+ fCandidatesTotal.push_back(cand);
fCandidates.push_back(cand);
continue;
}
@@ -1203,61 +1383,279 @@ void CbmAnaDielectronTask::FillCandidates()
IsElectron(iGTrack, cand.fMomentum.Mag(), &cand);
fCandidates.push_back(cand);
+ fCandidatesTotal.push_back(cand);
+
if (!cand.fIsElectron && cand.fChi2Prim < fCuts.fChiPrimCut) fTTCandidates.push_back(cand);
} // global tracks
cout << "fCandidates.size() = " << fCandidates.size() << endl;
cout << "fTTCandidates.size() = " << fTTCandidates.size() << endl;
- AssignMcToCandidates();
+ AssignMcToCandidates(fCandidates);
+ AssignMcToCandidates(fCandidatesTotal);
AssignMcToTopologyCandidates(fTTCandidates);
}
-void CbmAnaDielectronTask::AssignMcToCandidates()
+void CbmAnaDielectronTask::CombinatorialPairs()
{
- int nCand = fCandidates.size();
+ int nofEvents = fh_event_number->GetEntries();
+ cout << "Number of Events = " << nofEvents << endl;
+ int nCand = fCandidatesTotal.size();
+ int nSigCand1 = 0;
+ int nSigCand2 = 0;
+ int nChiPrim = 0;
+ int nIsEl = 0;
+ int nGammaCut = 0;
+ int nStCut = 0;
+ int nRtCut = 0;
+ int nTtCut = 0;
+ int nPtCut = 0;
+
+ for (int iCand1 = 0; iCand1 < nCand - 1; iCand1++) { // Before: iCand1 < nCand - 2
+ int nEvent1 = fCandidatesTotal[iCand1].fEventNumber;
+ int charge1 = fCandidatesTotal[iCand1].fCharge;
+ Bool_t isSignal1 = (fCandidatesTotal[iCand1].fIsMcSignalElectron);
+ if (isSignal1) nSigCand1++;
+
+ for (int iCand2 = iCand1 + 1; iCand2 < nCand; iCand2++) {
+ int nEvent2 = fCandidatesTotal[iCand2].fEventNumber;
+ int charge2 = fCandidatesTotal[iCand2].fCharge;
+ Bool_t isSignal2 = (fCandidatesTotal[iCand2].fIsMcSignalElectron);
+ if (isSignal2) nSigCand2++;
+ if (isSignal1 xor isSignal2)
+ continue; // since in FillPairHists() this case is not considered as well; to have same behaviour here as there
+ double weight = (isSignal1 || isSignal2) ? fWeight : 1.;
+ //double weight = 1;
+
+ CbmLmvmKinematicParams pRec =
+ CbmLmvmKinematicParams::KinematicParamsWithCandidates(&fCandidatesTotal[iCand1], &fCandidatesTotal[iCand2]);
+
+ Bool_t isChiPrimary = (fCandidatesTotal[iCand1].fChi2Prim < fCuts.fChiPrimCut
+ && fCandidatesTotal[iCand2].fChi2Prim < fCuts.fChiPrimCut);
+ Bool_t isEl = (fCandidatesTotal[iCand1].fIsElectron && fCandidatesTotal[iCand2].fIsElectron);
+ Bool_t isGammaCut = (!fCandidatesTotal[iCand1].fIsGamma && !fCandidatesTotal[iCand2].fIsGamma);
+ Bool_t isStCut = (fCandidatesTotal[iCand1].fIsStCutElectron && fCandidatesTotal[iCand2].fIsStCutElectron);
+ Bool_t isRtCut = (fCandidatesTotal[iCand1].fIsRtCutElectron && fCandidatesTotal[iCand2].fIsRtCutElectron);
+ Bool_t isTtCut = (fCandidatesTotal[iCand1].fIsTtCutElectron && fCandidatesTotal[iCand2].fIsTtCutElectron);
+ Bool_t isPtCut = (fCandidatesTotal[iCand1].fMomentum.Perp() > fCuts.fPtCut
+ && fCandidatesTotal[iCand2].fMomentum.Perp() > fCuts.fPtCut);
+ Bool_t isMvd1Cut = (fCandidatesTotal[iCand1].fIsMvd1CutElectron && fCandidatesTotal[iCand2].fIsMvd1CutElectron);
+ Bool_t isMvd2Cut = (fCandidatesTotal[iCand1].fIsMvd2CutElectron && fCandidatesTotal[iCand2].fIsMvd2CutElectron);
+
+ if (isChiPrimary) nChiPrim++;
+ if (isChiPrimary && isEl) nIsEl++;
+ if (isChiPrimary && isEl && isGammaCut) nGammaCut++;
+ if (isChiPrimary && isEl && isGammaCut && isStCut) nStCut++;
+ if (isChiPrimary && isEl && isGammaCut && isStCut && isRtCut) nRtCut++;
+ if (isChiPrimary && isEl && isGammaCut && isStCut && isRtCut && isTtCut) nTtCut++;
+ if (isChiPrimary && isEl && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut) nPtCut++;
+
+ if (!fUseMvd) isMvd1Cut = true;
+ if (!fUseMvd) isMvd2Cut = true;
+
+ // step: reco
+ if (charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kReco]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kReco]->Fill(pRec.fMinv, weight);
+ }
+ if (charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kReco]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kReco]->Fill(pRec.fMinv, weight);
+ }
+ if (charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kReco]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kReco]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: chi2prim
+ if (isChiPrimary && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kChi2Prim]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kChi2Prim]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kChi2Prim]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kChi2Prim]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kChi2Prim]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kChi2Prim]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: elID
+ if (isChiPrimary && isEl && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kElId]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kElId]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kElId]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kElId]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kElId]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kElId]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: gamma cut
+ if (isChiPrimary && isEl && isGammaCut && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kGammaCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kGammaCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kGammaCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kGammaCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kGammaCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kGammaCut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: MVD1 cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kMvd1Cut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: MVD2 cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kMvd2Cut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: ST cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kStCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kStCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kStCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kStCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kStCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kStCut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: RT cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kRtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kRtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && charge1 < 0
+ && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kRtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kRtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && charge1 > 0
+ && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kRtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kRtCut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: TT cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut
+ && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kTtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kTtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut && charge1 < 0
+ && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kTtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kTtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut && charge1 > 0
+ && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kTtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kTtCut]->Fill(pRec.fMinv, weight);
+ }
+
+ // step: Pt cut
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge1 * charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPM_minv_sameEvent[kPtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPM_minv_mixedEvents[kPtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge1 < 0 && charge2 < 0) {
+ if (nEvent1 == nEvent2) fh_combPairsMM_minv_sameEvent[kPtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsMM_minv_mixedEvents[kPtCut]->Fill(pRec.fMinv, weight);
+ }
+ if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut && isRtCut && isTtCut && isPtCut
+ && charge1 > 0 && charge2 > 0) {
+ if (nEvent1 == nEvent2) fh_combPairsPP_minv_sameEvent[kPtCut]->Fill(pRec.fMinv, weight);
+ if (nEvent1 != nEvent2) fh_combPairsPP_minv_mixedEvents[kPtCut]->Fill(pRec.fMinv, weight);
+ }
+
+ } // FOR cand2
+ } // FOR cand1
+
+ cout << "fCandidatesTotal.size() = " << fCandidatesTotal.size() << endl;
+ cout << "nSigCand1 = " << nSigCand1 << ", nSigCand2 = " << nSigCand2 << endl;
+ cout << "nChiPrim = " << nChiPrim << ", nIsEl = " << nIsEl << ", nGammaCut = " << nGammaCut
+ << ", nSt/Rt/TtCut = " << nStCut << "/" << nRtCut << "/" << nTtCut << ", nPtCut = " << nPtCut << endl;
+}
+
+void CbmAnaDielectronTask::AssignMcToCandidates(vector<CbmLmvmCandidate>& candVector)
+{
+ int nCand = candVector.size();
for (int i = 0; i < nCand; i++) {
- fCandidates[i].ResetMcParams();
+ if (candVector[i].fEventNumber != fEventNumber) continue;
+ candVector[i].ResetMcParams();
//STS
//MCTrackId of the candidate is defined by STS track
- int stsInd = fCandidates[i].fStsInd;
+ int stsInd = candVector[i].fStsInd;
CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
if (stsMatch == NULL) continue;
if (stsMatch->GetNofLinks() == 0) continue;
- fCandidates[i].fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
- if (fCandidates[i].fStsMcTrackId < 0) continue;
- CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMCTracks->At(fCandidates[i].fStsMcTrackId);
+ candVector[i].fStsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
+ if (candVector[i].fStsMcTrackId < 0) continue;
+ CbmMCTrack* mcTrack1 = (CbmMCTrack*) fMCTracks->At(candVector[i].fStsMcTrackId);
if (mcTrack1 == NULL) continue;
- int pdg = TMath::Abs(mcTrack1->GetPdgCode());
- int motherId = mcTrack1->GetMotherId();
- fCandidates[i].fMcMotherId = motherId;
- fCandidates[i].fMcPdg = pdg;
+ int pdg = TMath::Abs(mcTrack1->GetPdgCode());
+ int motherId = mcTrack1->GetMotherId();
+ candVector[i].fMcMotherId = motherId;
+ candVector[i].fMcPdg = pdg;
if (pdg == 211 && fPionMisidLevel >= 0.) continue;
// RICH
- int richInd = fCandidates[i].fRichInd;
+ int richInd = candVector[i].fRichInd;
CbmTrackMatchNew* richMatch = (CbmTrackMatchNew*) fRichRingMatches->At(richInd);
if (richMatch == NULL) continue;
- fCandidates[i].fRichMcTrackId = richMatch->GetMatchedLink().GetIndex();
+ candVector[i].fRichMcTrackId = richMatch->GetMatchedLink().GetIndex();
- fCandidates[i].fIsMcSignalElectron = CbmLmvmUtils::IsMcSignalElectron(mcTrack1);
- fCandidates[i].fIsMcPi0Electron = CbmLmvmUtils::IsMcPi0Electron(mcTrack1, fMCTracks);
- fCandidates[i].fIsMcGammaElectron = CbmLmvmUtils::IsMcGammaElectron(mcTrack1, fMCTracks);
- fCandidates[i].fIsMcEtaElectron = CbmLmvmUtils::IsMcEtaElectron(mcTrack1, fMCTracks);
+ candVector[i].fIsMcSignalElectron = CbmLmvmUtils::IsMcSignalElectron(mcTrack1);
+ candVector[i].fIsMcPi0Electron = CbmLmvmUtils::IsMcPi0Electron(mcTrack1, fMCTracks);
+ candVector[i].fIsMcGammaElectron = CbmLmvmUtils::IsMcGammaElectron(mcTrack1, fMCTracks);
+ candVector[i].fIsMcEtaElectron = CbmLmvmUtils::IsMcEtaElectron(mcTrack1, fMCTracks);
// TRD
// CbmTrdTrack* trdTrack = NULL;
if (fUseTrd == true) {
- int trdInd = fCandidates[i].fTrdInd;
+ int trdInd = candVector[i].fTrdInd;
CbmTrackMatchNew* trdMatch = (CbmTrackMatchNew*) fTrdTrackMatches->At(trdInd);
if (trdMatch == NULL) continue;
- fCandidates[i].fTrdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
+ candVector[i].fTrdMcTrackId = trdMatch->GetMatchedLink().GetIndex();
}
// ToF
- int tofInd = fCandidates[i].fTofInd;
+ int tofInd = candVector[i].fTofInd;
if (tofInd < 0) continue;
CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd);
if (tofHit == NULL) continue;
@@ -1267,8 +1665,8 @@ void CbmAnaDielectronTask::AssignMcToCandidates()
if (tofPointIndex < 0) continue;
FairMCPoint* tofPoint = (FairMCPoint*) fTofPoints->At(tofPointIndex);
if (tofPoint == NULL) continue;
- fCandidates[i].fTofMcTrackId = tofPoint->GetTrackID();
- } // candidates
+ candVector[i].fTofMcTrackId = tofPoint->GetTrackID();
+ }
}
void CbmAnaDielectronTask::AssignMcToTopologyCandidates(vector<CbmLmvmCandidate>& cutCandidates)
@@ -1445,10 +1843,12 @@ void CbmAnaDielectronTask::FillPairHists(CbmLmvmCandidate* candP, CbmLmvmCandida
Bool_t isBG = (!isEta) && (!isGamma) && (!isPi0) && (!(candP->fIsMcSignalElectron || candM->fIsMcSignalElectron));
Bool_t isMismatch = (IsMismatch(candP) || IsMismatch(candM));
- if (isSignal) fh_signal_pty[step]->Fill(parMc->fRapidity, parMc->fPt, fWeight);
- if (isSignal) fh_signal_mom[step]->Fill(parMc->fMomentumMag, fWeight);
- if (isSignal) fh_signal_minv[step]->Fill(parRec->fMinv, fWeight);
- if (isSignal) fh_signal_minv_pt[step]->Fill(parRec->fMinv, parMc->fPt, fWeight);
+ double weight = fWeight;
+
+ if (isSignal) fh_signal_pty[step]->Fill(parMc->fRapidity, parMc->fPt, weight);
+ if (isSignal) fh_signal_mom[step]->Fill(parMc->fMomentumMag, weight);
+ if (isSignal) fh_signal_minv[step]->Fill(parRec->fMinv, weight);
+ if (isSignal) fh_signal_minv_pt[step]->Fill(parRec->fMinv, parMc->fPt, weight);
if (isBG) fh_bg_minv[step]->Fill(parRec->fMinv);
PairSource(candP, candM, step, parRec);
if (isPi0) fh_pi0_minv[step]->Fill(parRec->fMinv);
@@ -1598,15 +1998,25 @@ void CbmAnaDielectronTask::CheckTopologyCut(const string& cutName, const vector<
{
vector<Double_t> angles, mom, candInd;
Int_t nCand = fCandidates.size();
+ Int_t nCandTot = fCandidatesTotal.size();
Int_t nCutCand = cutCandidates.size();
- for (Int_t iP = 0; iP < nCand; iP++) {
+ int nJump = 0;
+ for (int i = 0; i < nCandTot; i++) {
+ if (fCandidatesTotal[i].fEventNumber != fEventNumber) nJump++;
+ if (fCandidatesTotal[i].fEventNumber == fEventNumber) {
+ cout << "nJump = " << nJump << endl;
+ break;
+ }
+ }
+
+ for (Int_t iP = 0; iP < nCand; iP++) {
if (fCandidates[iP].fChi2Prim < fCuts.fChiPrimCut && fCandidates[iP].fIsElectron) {
angles.clear();
mom.clear();
candInd.clear();
for (Int_t iM = 0; iM < nCutCand; iM++) {
- // different charges, charge Im != charge iP
+ // different charges, charge iM != charge iP
if (cutCandidates[iM].fCharge != fCandidates[iP].fCharge) {
CbmLmvmKinematicParams pRec =
CbmLmvmKinematicParams::KinematicParamsWithCandidates(&fCandidates[iP], &cutCandidates[iM]);
@@ -1625,17 +2035,35 @@ void CbmAnaDielectronTask::CheckTopologyCut(const string& cutName, const vector<
}
}
if (minInd == -1) {
- if (cutName == "TT") fCandidates[iP].fIsTtCutElectron = true;
- if (cutName == "ST") fCandidates[iP].fIsStCutElectron = true;
- if (cutName == "RT") fCandidates[iP].fIsRtCutElectron = true;
+ if (cutName == "TT") {
+ fCandidates[iP].fIsTtCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsTtCutElectron = true;
+ }
+ if (cutName == "ST") {
+ fCandidates[iP].fIsStCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsStCutElectron = true;
+ }
+ if (cutName == "RT") {
+ fCandidates[iP].fIsRtCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsRtCutElectron = true;
+ }
continue;
}
Double_t sqrt_mom = TMath::Sqrt(fCandidates[iP].fMomentum.Mag() * mom[minInd]);
Double_t val = -1. * (angleCut / ppCut) * sqrt_mom + angleCut;
if (!(sqrt_mom < ppCut && val > minAng)) {
- if (cutName == "TT") fCandidates[iP].fIsTtCutElectron = true;
- if (cutName == "ST") fCandidates[iP].fIsStCutElectron = true;
- if (cutName == "RT") fCandidates[iP].fIsRtCutElectron = true;
+ if (cutName == "TT") {
+ fCandidates[iP].fIsTtCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsTtCutElectron = true;
+ }
+ if (cutName == "ST") {
+ fCandidates[iP].fIsStCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsStCutElectron = true;
+ }
+ if (cutName == "RT") {
+ fCandidates[iP].fIsRtCutElectron = true;
+ fCandidatesTotal[iP + nJump].fIsRtCutElectron = true;
+ }
}
int stsInd = cutCandidates[candInd[minInd]].fStsInd;
@@ -1671,7 +2099,7 @@ void CbmAnaDielectronTask::CheckTopologyCut(const string& cutName, const vector<
}
} //if electron
} //iP
- /*
+ /*
cout << "Opening angles between cand and cutCand: " << endl;
for(int i = 0; i < angles.size(); i++){
cout << angles[i] << ", ";
@@ -1793,7 +2221,7 @@ void CbmAnaDielectronTask::IsElectron(Int_t globalTrackIndex, Double_t momentum,
}
else {
// PID using MC information, a certain pi supression level can be set
- if (cand->fStsMcTrackId < 0 || cand->fStsMcTrackId >= fMCTracks->GetEntriesFast()) { cand->fIsElectron = false; }
+ if (cand->fStsMcTrackId < 0 || cand->fStsMcTrackId >= fMCTracks->GetEntries()) { cand->fIsElectron = false; }
else {
CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(cand->fStsMcTrackId);
Int_t pdg = mcTrack->GetPdgCode();
@@ -1809,7 +2237,6 @@ void CbmAnaDielectronTask::IsElectron(Int_t globalTrackIndex, Double_t momentum,
}
}
-
void CbmAnaDielectronTask::DifferenceSignalAndBg()
{
//ID cuts
@@ -2115,6 +2542,7 @@ void CbmAnaDielectronTask::MvdCutMcDistance()
void CbmAnaDielectronTask::Finish()
{
+ CombinatorialPairs();
TDirectory* oldir = gDirectory;
TFile* outFile = FairRootManager::Instance()->GetOutFile();
if (outFile != NULL) {
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
old mode 100644
new mode 100755
index 373af134226db07ff52b867db259c45eccfc8c7e..d46ec889f487e405cbc85a29770979ddbc788635
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
@@ -6,7 +6,7 @@
* @author Elena Lebedeva <e.lebedeva@gsi.de>
* @since 2010
* @version 3.0
- **/
+**/
#ifndef CBM_ANA_DIELECTRON_TASK_H
#define CBM_ANA_DIELECTRON_TASK_H
@@ -140,7 +140,7 @@ public:
CbmLmvmKinematicParams* parRec);
/*
- * \brief Fills minv, pty, mom histograms for specified analysis step.
+ * \brief Fills minv, pty, mom histograms for specified analysis step.
* \param[in] candP Positive candidate.
* \param[in] candM Negative candidate.
* \param[in] parMc MC kinematic parameters.
@@ -168,8 +168,9 @@ public:
void FillCandidates();
-
- void AssignMcToCandidates();
+ void AssignMcToCandidates(std::vector<CbmLmvmCandidate>& candVector);
+ //void AssignMcToCandidates();
+ //void AssignMcToCandidatesTotal();
void AssignMcToTopologyCandidates(std::vector<CbmLmvmCandidate>& cutCandidates);
@@ -184,7 +185,7 @@ public:
void CheckGammaConvAndPi0();
- void FillNofChargedParticlesPerEvent();
+ void FillNofChargedParticles();
/*
* \brief
@@ -201,6 +202,12 @@ public:
const std::vector<TH2D*>& hcut, const std::vector<TH2D*>& hcutPion,
const std::vector<TH2D*>& hcutTruepair, Double_t angleCut, Double_t ppCut);
+ /*
+ * \brief Set cut values and fill histograms for topology cut
+ * \param[in] cutName ST or TT
+ */
+ // void CheckTopologyCutTotal(const std::string& cutName, const std::vector<CbmLmvmCandidate>& cutCandidates, Double_t angleCut, Double_t ppCut);
+
void CalculateNofTopologyPairs(TH1D* h_nof_pairs, const std::string& source);
void MvdCutMcDistance();
@@ -208,6 +215,8 @@ public:
// Likelihood vs Momentum
void FillMomLikeHist();
+ void CombinatorialPairs();
+
virtual void Finish();
void FillElPiMomHist();
@@ -257,6 +266,7 @@ private:
Bool_t fUseTof;
std::vector<CbmLmvmCandidate> fCandidates;
+ std::vector<CbmLmvmCandidate> fCandidatesTotal;
std::vector<CbmLmvmCandidate> fSTCandidates; // STCut Segmented tracks, reconstructed only in STS
std::vector<CbmLmvmCandidate>
fTTCandidates; // TTCut Reconstructed tracks, reconstructed in all detectors but not identified as electrons
@@ -268,6 +278,8 @@ private:
Double_t fPionMisidLevel; // For the ideal particle identification cases, set to -1 for real PID
TRandom3* fRandom3;
+ Int_t fEventNumber; // number of current event
+
//Bool_t fUseMcMomentum;
CbmLmvmCuts fCuts; // electorn identification and analisys cuts
@@ -279,8 +291,8 @@ private:
TH2D* fh_mc_signal_mom_angle; // angle vs. sqrt(mom1*mom2) with MCTracks
- TH1D* fh_nof_charged_particles; //charged UrQMD particles
- TH1D* fh_nof_charged_particles_acc; //accepted charged UrQMD particles
+ TH1D* fh_nof_charged_particles; // charged UrQMD particles
+ TH1D* fh_nof_charged_particles_acc; // accepted charged UrQMD particles
TH1D* fh_mc_mother_pdg; //mother pdg code for e-/e+
TH1D* fh_acc_mother_pdg; //mother pdg code for accepted e-/e+
@@ -302,8 +314,20 @@ private:
// [5]-gamma cut, [6]-mvd1cut, [7]-mvd2cut, [8]-stcut, [9]-ttcut, [10]-ptcut.
//Use AnalysisSteps enumeration for access.
//MC and ACC histograms are not filled sometimes.
- std::vector<TH1D*> fh_signal_minv; // Invariant mass for Signal
- std::vector<TH1D*> fh_bg_minv; // Invariant mass for BG
+ std::vector<TH1D*> fh_signal_minv; // Invariant mass for Signal
+ std::vector<TH1D*> fh_bg_minv; // Invariant mass for BG
+
+ std::vector<TH1D*> fh_combPairsPM_minv_sameEvent; // Invariant mass for comb. Pairs of e+/e-
+ std::vector<TH1D*> fh_combPairsPP_minv_sameEvent; // Invariant mass for comb. Pairs of e+/e+
+ std::vector<TH1D*> fh_combPairsMM_minv_sameEvent; // Invariant mass for comb. Pairs of e-/e-
+ std::vector<TH1D*> fh_combPairsPM_minv_mixedEvents; // Invariant mass for comb. Pairs of e+/e-
+ std::vector<TH1D*> fh_combPairsPP_minv_mixedEvents; // Invariant mass for comb. Pairs of e+/e+
+ std::vector<TH1D*> fh_combPairsMM_minv_mixedEvents; // Invariant mass for comb. Pairs of e-/e-
+ std::vector<TH1D*> fh_nof_plutoElectrons; // nof electrons / positrons vs. momentum
+ std::vector<TH1D*> fh_nof_plutoPositrons; // nof electrons / positrons vs. momentum
+ std::vector<TH1D*> fh_nof_urqmdElectrons; // nof electrons / positrons vs. momentum
+ std::vector<TH1D*> fh_nof_urqmdPositrons; // nof electrons / positrons vs. momentum
+
std::vector<TH1D*> fh_pi0_minv; // Invariant mass for Pi0
std::vector<TH1D*> fh_eta_minv; // Invariant mass for Eta
std::vector<TH1D*> fh_gamma_minv; // Invariant mass for Eta
@@ -374,6 +398,7 @@ private:
//store event number
TH1D* fh_event_number;
+ TH1D* fh_event_number_mixed; // number of involved mixed events
//nof signal and bg tracks after each cut
TH1D* fh_nof_bg_tracks;
@@ -442,7 +467,6 @@ public:
void SetUseTof(Bool_t use) { fUseTof = use; };
void SetWeight(Double_t weight) { fWeight = weight; };
void SetEnergyAndPlutoParticle(const string& energy, const string& particle);
-
void SetPionMisidLevel(Double_t level) { fPionMisidLevel = level; }
// void SetMomentumCut(Double_t mom) {fMomentumCut = mom;}
};
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
old mode 100644
new mode 100755
index d8f8a02494c53274f4965b527a2a525541e90886..c3d3b24ad4a82f406a6e49a28b2b14ded0681659
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
@@ -70,7 +70,6 @@ void CbmAnaDielectronTaskDraw::DrawHistFromFile(const string& fileName, const st
TFile* file = new TFile(fileName.c_str());
fHM->ReadFromFile(file);
fNofEvents = (Int_t) H1("fh_event_number")->GetEntries();
- //fNofEvents = 10000;
cout << "File name = " << fileName << endl;
cout << "Number of events = " << fNofEvents << endl;
@@ -100,6 +99,7 @@ void CbmAnaDielectronTaskDraw::DrawHistFromFile(const string& fileName, const st
DrawElPiMomHis();
DrawPmtXY();
DrawMomLikeHist();
+ DrawSingleParticleYield();
SaveCanvasToImage();
/// Restore old global file and folder pointer to avoid messing with FairRoot
@@ -118,12 +118,144 @@ void CbmAnaDielectronTaskDraw::DrawMomLikeHist()
DrawH2(H2("fh_mom_likelihood_Pi"));
}
+void CbmAnaDielectronTaskDraw::DrawSingleParticleYield()
+{
+ TCanvas* c2 = fHM->CreateCanvas("fh_mom_ElPos_pluto", "fh_mom_ElPos_pluto", 1000, 500);
+ c2->Divide(2, 1);
+ c2->cd(1);
+ Draw1DHistoForEachAnalysisStep("fh_nof_plutoElectrons", true);
+ c2->cd(2);
+ Draw1DHistoForEachAnalysisStep("fh_nof_plutoPositrons", true);
+
+ TCanvas* c3 = fHM->CreateCanvas("fh_mom_ElPos_urqmd", "fh_mom_ElPos_urqmd", 1000, 500);
+ c3->Divide(2, 1);
+ c3->cd(1);
+ Draw1DHistoForEachAnalysisStep("fh_nof_urqmdElectrons", true);
+ c3->cd(2);
+ Draw1DHistoForEachAnalysisStep("fh_nof_urqmdPositrons", true);
+
+ // draw electron and positron yield vs. momentum
+ TH1D* nPlutoElMc = (TH1D*) H1("fh_nof_plutoElectrons_mc")->Clone();
+ TH1D* nPlutoPosMc = (TH1D*) H1("fh_nof_plutoPositrons_mc")->Clone();
+ TH1D* nPlutoElAcc = (TH1D*) H1("fh_nof_plutoElectrons_acc")->Clone();
+ TH1D* nPlutoPosAcc = (TH1D*) H1("fh_nof_plutoPositrons_acc")->Clone();
+ TH1D* nPlutoElReco = (TH1D*) H1("fh_nof_plutoElectrons_reco")->Clone();
+ TH1D* nPlutoPosReco = (TH1D*) H1("fh_nof_plutoPositrons_reco")->Clone();
+ TH1D* nPlutoElElid = (TH1D*) H1("fh_nof_plutoElectrons_elid")->Clone();
+ TH1D* nPlutoPosElid = (TH1D*) H1("fh_nof_plutoPositrons_elid")->Clone();
+ TH1D* nPlutoElTt = (TH1D*) H1("fh_nof_plutoElectrons_ttcut")->Clone();
+ TH1D* nPlutoPosTt = (TH1D*) H1("fh_nof_plutoPositrons_ttcut")->Clone();
+ TH1D* nPlutoElPt = (TH1D*) H1("fh_nof_plutoElectrons_ptcut")->Clone();
+ TH1D* nPlutoPosPt = (TH1D*) H1("fh_nof_plutoPositrons_ptcut")->Clone();
+ TH1D* nUrqmdElMc = (TH1D*) H1("fh_nof_urqmdElectrons_mc")->Clone();
+ TH1D* nUrqmdPosMc = (TH1D*) H1("fh_nof_urqmdPositrons_mc")->Clone();
+ TH1D* nUrqmdElAcc = (TH1D*) H1("fh_nof_urqmdElectrons_acc")->Clone();
+ TH1D* nUrqmdPosAcc = (TH1D*) H1("fh_nof_urqmdPositrons_acc")->Clone();
+ TH1D* nUrqmdElReco = (TH1D*) H1("fh_nof_urqmdElectrons_reco")->Clone();
+ TH1D* nUrqmdPosReco = (TH1D*) H1("fh_nof_urqmdPositrons_reco")->Clone();
+ TH1D* nUrqmdElElid = (TH1D*) H1("fh_nof_urqmdElectrons_elid")->Clone();
+ TH1D* nUrqmdPosElid = (TH1D*) H1("fh_nof_urqmdPositrons_elid")->Clone();
+ TH1D* nUrqmdElTt = (TH1D*) H1("fh_nof_urqmdElectrons_ttcut")->Clone();
+ TH1D* nUrqmdPosTt = (TH1D*) H1("fh_nof_urqmdPositrons_ttcut")->Clone();
+ TH1D* nUrqmdElPt = (TH1D*) H1("fh_nof_urqmdElectrons_ptcut")->Clone();
+ TH1D* nUrqmdPosPt = (TH1D*) H1("fh_nof_urqmdPositrons_ptcut")->Clone();
+
+ double min1 = 5e-8;
+ double max1 = 50;
+ nPlutoElMc->SetMinimum(min1);
+ nPlutoElMc->SetMaximum(max1);
+ nPlutoPosMc->SetMinimum(min1);
+ nPlutoPosMc->SetMaximum(max1);
+ nPlutoElAcc->SetMinimum(min1);
+ nPlutoElAcc->SetMaximum(max1);
+ nPlutoPosAcc->SetMinimum(min1);
+ nPlutoPosAcc->SetMaximum(max1);
+ nPlutoElReco->SetMinimum(min1);
+ nPlutoElReco->SetMaximum(max1);
+ nPlutoPosReco->SetMinimum(min1);
+ nPlutoPosReco->SetMaximum(max1);
+ nPlutoElElid->SetMinimum(min1);
+ nPlutoElElid->SetMaximum(max1);
+ nPlutoPosElid->SetMinimum(min1);
+ nPlutoPosElid->SetMaximum(max1);
+ nPlutoElTt->SetMinimum(min1);
+ nPlutoElTt->SetMaximum(max1);
+ nPlutoPosTt->SetMinimum(min1);
+ nPlutoPosTt->SetMaximum(max1);
+ nPlutoElPt->SetMinimum(min1);
+ nPlutoElPt->SetMaximum(max1);
+ nPlutoPosPt->SetMinimum(min1);
+ nPlutoPosPt->SetMaximum(max1);
+ nUrqmdElMc->SetMinimum(min1);
+ nUrqmdElMc->SetMaximum(max1);
+ nUrqmdPosMc->SetMinimum(min1);
+ nUrqmdPosMc->SetMaximum(max1);
+ nUrqmdElAcc->SetMinimum(min1);
+ nUrqmdElAcc->SetMaximum(max1);
+ nUrqmdPosAcc->SetMinimum(min1);
+ nUrqmdPosAcc->SetMaximum(max1);
+ nUrqmdElReco->SetMinimum(min1);
+ nUrqmdElReco->SetMaximum(max1);
+ nUrqmdPosReco->SetMinimum(min1);
+ nUrqmdPosReco->SetMaximum(max1);
+ nUrqmdElElid->SetMinimum(min1);
+ nUrqmdElElid->SetMaximum(max1);
+ nUrqmdPosElid->SetMinimum(min1);
+ nUrqmdPosElid->SetMaximum(max1);
+ nUrqmdElTt->SetMinimum(min1);
+ nUrqmdElTt->SetMaximum(max1);
+ nUrqmdPosTt->SetMinimum(min1);
+ nUrqmdPosTt->SetMaximum(max1);
+ nUrqmdElPt->SetMinimum(min1);
+ nUrqmdElPt->SetMaximum(max1);
+ nUrqmdPosPt->SetMinimum(min1);
+ nUrqmdPosPt->SetMaximum(max1);
+
+ fHM->CreateCanvas("fh_nof_plutoElPos", "fh_nof_plutoElPos", 800, 800);
+ DrawH1({nPlutoElMc, nPlutoPosMc, nPlutoElAcc, nPlutoPosAcc, nPlutoElReco, nPlutoPosReco}, {"", "", "", "", "", ""},
+ kLinear, kLog, false, 0, 0, 0, 0, "Hist p");
+
+ TLegend* legendNofPP = new TLegend(0.65, 0.6, 0.88, 0.93);
+ legendNofPP->SetFillColor(kWhite);
+ legendNofPP->AddEntry(nPlutoElMc, "electrons kMc");
+ legendNofPP->AddEntry(nPlutoPosMc, "positrons kMc");
+ legendNofPP->AddEntry(nPlutoElAcc, "electrons kAcc");
+ legendNofPP->AddEntry(nPlutoPosAcc, "positrons kAcc");
+ legendNofPP->AddEntry(nPlutoElReco, "electrons kReco");
+ legendNofPP->AddEntry(nPlutoPosReco, "positrons kReco");
+ legendNofPP->Draw();
+
+ fHM->CreateCanvas("fh_nof_urqmdElPos", "fh_nof_urqmdElPos", 800, 800);
+ DrawH1({nUrqmdElMc, nUrqmdPosMc, nUrqmdElAcc, nUrqmdPosAcc, nUrqmdElReco, nUrqmdPosReco}, {"", "", "", "", "", ""},
+ kLinear, kLog, false, 0, 0, 0, 0, "Hist p");
+
+ TLegend* legendNofUP = new TLegend(0.65, 0.6, 0.88, 0.93);
+ legendNofUP->SetFillColor(kWhite);
+ legendNofUP->AddEntry(nUrqmdElMc, "electrons kMc");
+ legendNofUP->AddEntry(nUrqmdPosMc, "positrons kMc");
+ legendNofUP->AddEntry(nUrqmdElAcc, "electrons kAcc");
+ legendNofUP->AddEntry(nUrqmdPosAcc, "positrons kAcc");
+ legendNofUP->AddEntry(nUrqmdElReco, "electrons kReco");
+ legendNofUP->AddEntry(nUrqmdPosReco, "positrons kReco");
+ legendNofUP->Draw();
+}
+
void CbmAnaDielectronTaskDraw::RebinMinvHist()
{
- int nRebin = 10;
+ int nRebin = 20;
for (int i = 0; i < CbmLmvmHist::fNofAnaSteps; i++) {
H1("fh_signal_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
H1("fh_bg_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsPM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsPP_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsMM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsPM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsPP_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_combPairsMM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
H1("fh_pi0_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
H1("fh_eta_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
H1("fh_bg_truematch_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(2 * nRebin);
@@ -131,6 +263,15 @@ void CbmAnaDielectronTaskDraw::RebinMinvHist()
H1("fh_bg_truematch_el_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(2 * nRebin);
H1("fh_bg_truematch_notel_minv_" + CbmLmvmHist::fAnaSteps[i])->Rebin(2 * nRebin);
+ H1("fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebin);
+ H1("fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebin);
+ H1("fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebin);
+ H1("fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebin);
+ H1("fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebin);
+
for (int iP = 0; iP < CbmLmvmHist::fNofBgPairSources; iP++) {
stringstream ss;
ss << "fh_source_bg_minv_" << iP << "_" << CbmLmvmHist::fAnaSteps[i];
@@ -722,6 +863,7 @@ void CbmAnaDielectronTaskDraw::DrawMinvForEachAnalysisStep()
c1->cd(2);
Draw1DHistoForEachAnalysisStep("fh_bg_minv", true);
+
TCanvas* c2 = fHM->CreateCanvas("lmvm_minv_for_each_analysis_step_pi0_eta",
"lmvm_minv_for_each_analysis_step_pi0_eta", 1200, 600);
c2->Divide(2, 1);
@@ -730,6 +872,26 @@ void CbmAnaDielectronTaskDraw::DrawMinvForEachAnalysisStep()
c2->cd(2);
Draw1DHistoForEachAnalysisStep("fh_eta_minv", true);
+ TCanvas* c3 = fHM->CreateCanvas("lmvm_minv_for_each_analysis_step_combinatorial_Pairs_same_Event",
+ "lmvm_minv_for_each_analysis_step_combinatorial_Pairs_same_event", 1800, 600);
+ c3->Divide(3, 1);
+ c3->cd(1);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsPM_minv_sameEvent", true);
+ c3->cd(2);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsPP_minv_sameEvent", true);
+ c3->cd(3);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsMM_minv_sameEvent", true);
+
+ TCanvas* c4 = fHM->CreateCanvas("lmvm_minv_for_each_analysis_step_combinatorial_Pairs_mixed_Events",
+ "lmvm_minv_for_each_analysis_step_combinatorial_Pairs_mixed_Events", 1800, 600);
+ c4->Divide(3, 1);
+ c4->cd(1);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsPM_minv_mixedEvents", true);
+ c4->cd(2);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsPP_minv_mixedEvents", true);
+ c4->cd(3);
+ Draw1DHistoForEachAnalysisStep("fh_combPairsMM_minv_mixedEvents", true);
+
fHM->CreateCanvas("lmvm_minv_for_each_analysis_step_gamma", "lmvm_minv_for_each_analysis_step_gamma", 600, 600);
// H1("fh_gamma_minv_mc")->GetXaxis()->SetRangeUser(0., 0.05);
Draw1DHistoForEachAnalysisStep("fh_gamma_minv", true);
@@ -1242,4 +1404,7 @@ void CbmAnaDielectronTaskDraw::DrawMvdAndStsHist()
}
-void CbmAnaDielectronTaskDraw::SaveCanvasToImage() { fHM->SaveCanvasToImage(fOutputDir, "png;eps"); }
+void CbmAnaDielectronTaskDraw::SaveCanvasToImage()
+{
+ fHM->SaveCanvasToImage(fOutputDir, "png"); // fHM->SaveCanvasToImage(fOutputDir, "png;eps");
+}
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.h b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.h
old mode 100644
new mode 100755
index da7f3bc446389beef6b05f175d35961bf9f15020..63214525e80490d3667ec3af943ab79f3248b2c9
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.h
@@ -46,7 +46,8 @@ public:
Bool_t drawSig = true);
private:
- Int_t fNofEvents;
+ Int_t fNofEvents; // number of events of current job
+
Bool_t fUseMvd; // do you want to draw histograms related to the MVD detector?
Bool_t fDrawSignificance; // do you want to draw significance histograms of 1D cuts?
@@ -237,6 +238,9 @@ private:
// Draw Likelihood vs Momentum
void DrawMomLikeHist();
+ // Draw yield of electrons and positrons vs. momentum
+ void DrawSingleParticleYield();
+
CbmAnaDielectronTaskDraw(const CbmAnaDielectronTaskDraw&);
CbmAnaDielectronTaskDraw& operator=(const CbmAnaDielectronTaskDraw&);
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
index 9cae31ca15621b4ad40e0e88d0550c459c2b032d..5b720d2571c439a8ee34059a16812755e6d29b0d 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
@@ -46,6 +46,8 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
fUseMvd = useMvd;
fDrawQgp = (fileNameQgp != "");
+ nRebin = 100;
+
//SetDefaultDrawStyle();
vector<string> fileNames = {fileNameInmed, fileNameQgp, fileNameOmega, fileNamePhi, fileNameOmegaDalitz};
@@ -60,7 +62,7 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
TFile* file = new TFile(fileNames[i].c_str());
fHM[i]->ReadFromFile(file);
Int_t nofEvents = (int) H1(i, "fh_event_number")->GetEntries();
- fHM[i]->ScaleByPattern(".*", 1. / nofEvents);
+ //fHM[i]->ScaleByPattern(".*", 1. / nofEvents); // TODO: keep this commented?
cout << "nofEvents = " << nofEvents << endl;
}
@@ -72,17 +74,49 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
fh_mean_eta_minv_pt.resize(CbmLmvmHist::fNofAnaSteps);
fh_mean_pi0_minv_pt.resize(CbmLmvmHist::fNofAnaSteps);
fh_mean_sbg_vs_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPM_sameEvent_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPP_sameEvent_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsMM_sameEvent_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPM_mixedEvents_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPP_mixedEvents_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsMM_mixedEvents_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPM_sameEvent_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPP_sameEvent_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsMM_sameEvent_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPM_mixedEvents_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsPP_mixedEvents_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combPairsMM_mixedEvents_minv_raw.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_raw_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_assemb_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_GeomMeanSame_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_GeomMeanMixed_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_k_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combSignalNpm_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combSignalNpm_assemb_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ //fh_mean_combSignalBCoc_minv.resize(CbmLmvmHist::fNofAnaSteps); // needed?
+ fh_mean_combSignalBCoc_assemb_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combBg_errProp_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combSignal_errProp_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_combSBg_vs_minv.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_nof_plutoElectrons.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_nof_plutoPositrons.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_nof_urqmdElectrons.resize(CbmLmvmHist::fNofAnaSteps);
+ fh_mean_nof_urqmdPositrons.resize(CbmLmvmHist::fNofAnaSteps);
FillMeanHist();
FillSumSignalsHist();
CalcCutEffRange(0.0, 0.2);
CalcCutEffRange(0.2, 0.6);
CalcCutEffRange(0.6, 1.2);
+ CalcCombBGHistos();
SBgRangeAll();
DrawSBgSignals();
DrawMinvAll();
+ DrawMinvCombSignalAndBg();
DrawMinvPtAll();
DrawSBgVsMinv();
+ //CompareSTSversions();
SaveHist();
SaveCanvasToImage();
@@ -106,6 +140,20 @@ TH1D* CbmAnaDielectronTaskDrawAll::GetCoctailMinv(CbmLmvmAnalysisSteps step)
TH1D* sPi0 = fh_mean_pi0_minv[step];
TH1D* sOmegaDalitz = (TH1D*) H1(kOmegaD, "fh_signal_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ Int_t nofEvents = 0;
+ for (int i = 0; i < fNofSignals; i++) {
+ if (!fDrawQgp && i == kQgp) continue;
+ Int_t nofEventsS = (int) H1(i, "fh_event_number")->GetEntries();
+ if (i == 0) sInmed->Scale(1. / nofEventsS);
+ if (i == 1) sQgp->Scale(1. / nofEventsS);
+ if (i == 2) sOmega->Scale(1. / nofEventsS);
+ if (i == 3) sPhi->Scale(1. / nofEventsS);
+ if (i == 4) sOmegaDalitz->Scale(1. / nofEventsS);
+
+ nofEvents += nofEventsS;
+ if (i == 4) cout << "GetCoctailMinv(...): nofEvents = " << nofEvents << endl;
+ }
+
TH1D* coctail = (TH1D*) sInmed->Clone();
if (fDrawQgp) coctail->Add(sQgp);
coctail->Add(sOmega);
@@ -114,6 +162,8 @@ TH1D* CbmAnaDielectronTaskDrawAll::GetCoctailMinv(CbmLmvmAnalysisSteps step)
coctail->Add(sPi0);
coctail->Add(sOmegaDalitz);
+ //coctail->Scale(1. / (fNofSignals + 2)); // '+2' because except signals there are two particles (eta and pi0)
+
return coctail;
}
@@ -146,10 +196,24 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
TH1D* sPi0 = (TH1D*) fh_mean_pi0_minv[step]->Clone();
TH1D* sOmegaDalitz = (TH1D*) H1(kOmegaD, "fh_signal_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ Int_t nofEvents = 0;
+ for (int i = 0; i < fNofSignals; i++) {
+ if (!fDrawQgp && i == kQgp) continue;
+ Int_t nofEventsS = (int) H1(i, "fh_event_number")->GetEntries();
+ if (i == 0) sInmed->Scale(1. / nofEventsS);
+ if (i == 1) sQgp->Scale(1. / nofEventsS);
+ if (i == 2) sOmega->Scale(1. / nofEventsS);
+ if (i == 3) sPhi->Scale(1. / nofEventsS);
+ if (i == 4) sOmegaDalitz->Scale(1. / nofEventsS);
+
+ nofEvents += nofEventsS;
+ cout << "DrawMinv(...): nofEvents = " << nofEvents << endl;
+ }
+
TH1D* coctail = GetCoctailMinv(step);
TH1D* sbg = (TH1D*) bg->Clone();
- sbg->Add(sInmed);
+ sbg->Add(sInmed); // Why here not added coctail instead (is same procedure)?
if (fDrawQgp) sbg->Add(sQgp);
sbg->Add(sOmega);
sbg->Add(sPhi);
@@ -157,8 +221,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
sbg->Add(sPi0);
sbg->Add(sOmegaDalitz);
-
- int nRebin = 20;
sbg->Rebin(nRebin);
coctail->Rebin(nRebin);
bg->Rebin(nRebin);
@@ -170,17 +232,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
if (fDrawQgp) sQgp->Rebin(nRebin);
sPhi->Rebin(nRebin);
- /*sbg->Scale(1./nRebin);
- coctail->Scale(1./nRebin);
- bg->Scale(1./nRebin);
- sPi0->Scale(1./nRebin);
- sEta->Scale(1./nRebin);
- sOmegaDalitz->Scale(1./nRebin);
- sOmega->Scale(1./nRebin);
- sInmed->Scale(1./nRebin);
- sQgp->Scale(1./nRebin);
- sPhi->Scale(1./nRebin);*/
-
double binWidth = sbg->GetBinWidth(1);
sbg->Scale(1. / binWidth);
coctail->Scale(1. / binWidth);
@@ -193,7 +244,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
if (fDrawQgp) sQgp->Scale(1. / binWidth);
sPhi->Scale(1. / binWidth);
-
sbg->SetMinimum(5e-8);
sbg->SetMaximum(2e-2);
sbg->GetXaxis()->SetRangeUser(0, 2.);
@@ -258,7 +308,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
sQgp->SetFillStyle(3444);
}
-
sOmega->SetFillColor(kOrange + 7);
sOmega->SetLineColor(kOrange + 4);
sOmega->SetLineStyle(0);
@@ -315,9 +364,9 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
if (fDrawQgp) legend->AddEntry(sQgp, "QGP radiation", "f");
legend->AddEntry(sEta, "#eta #rightarrow #gammae^{+}e^{-}", "f");
legend->AddEntry(sPi0, "#pi^{0} #rightarrow #gammae^{+}e^{-}", "f");
- legend->AddEntry(coctail, "Coctail", "f");
+ legend->AddEntry(coctail, "Cocktail", "f");
legend->AddEntry(bg, "Background", "f");
- legend->AddEntry(sbg, "Coctail+BG", "f");
+ legend->AddEntry(sbg, "Cocktail+BG", "f");
legend->Draw();
}
else {
@@ -331,18 +380,1006 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
if (fDrawQgp) legend->AddEntry(sQgp, "QGP radiation", "f");
legend->AddEntry(sEta, "#eta #rightarrow #gammae^{+}e^{-}", "f");
legend->AddEntry(sPi0, "#pi^{0} #rightarrow #gammae^{+}e^{-}", "f");
- legend->AddEntry(coctail, "Coctail", "f");
+ legend->AddEntry(coctail, "Cocktail", "f");
legend->Draw();
}
+ gPad->SetLogy(true);
+}
+void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
+{
+ /******************************************************************************
+ * MIND: *
+ * 'B' (capital) stands for same event, 'b' (not-capital) for mixed events! *
+ * *
+ * Order of draw functions: *
+ * 1) Yields of single e- and e+ per event *
+ * 2) e+e- | e+e+ | e-e- pairs from same and mixed events *
+ * 3) geometric mean (conv. and assembled) *
+ * 4) B_comb with MC BG and B+- *
+ * 5) k factor *
+ * 6) comb. signal *
+ ******************************************************************************/
+
+ double yMin = 1e-7;
+ double yMax = 5e-2;
+ string yTitle = "dN/dM_{ee} [GeV/c^{2}]^{-1}";
+ string xTitle = "M_{ee} [GeV/c^2]";
+
+ bool setMinMax = false;
+
+ /* 1) draw number of e+ and e- per event */
+ // vs. momentum
+ TH1D* nPlutoElMc = (TH1D*) fh_mean_nof_plutoElectrons[kMc]->Clone();
+ TH1D* nPlutoPosMc = (TH1D*) fh_mean_nof_plutoPositrons[kMc]->Clone();
+ TH1D* nPlutoElAcc = (TH1D*) fh_mean_nof_plutoElectrons[kAcc]->Clone();
+ TH1D* nPlutoPosAcc = (TH1D*) fh_mean_nof_plutoPositrons[kAcc]->Clone();
+ TH1D* nPlutoElReco = (TH1D*) fh_mean_nof_plutoElectrons[kReco]->Clone();
+ TH1D* nPlutoPosReco = (TH1D*) fh_mean_nof_plutoPositrons[kReco]->Clone();
+ TH1D* nPlutoElElid = (TH1D*) fh_mean_nof_plutoElectrons[kElId]->Clone();
+ TH1D* nPlutoPosElid = (TH1D*) fh_mean_nof_plutoPositrons[kElId]->Clone();
+ TH1D* nPlutoElTt = (TH1D*) fh_mean_nof_plutoElectrons[kTtCut]->Clone();
+ TH1D* nPlutoPosTt = (TH1D*) fh_mean_nof_plutoPositrons[kTtCut]->Clone();
+ TH1D* nPlutoElPt = (TH1D*) fh_mean_nof_plutoElectrons[kPtCut]->Clone();
+ TH1D* nPlutoPosPt = (TH1D*) fh_mean_nof_plutoPositrons[kPtCut]->Clone();
+ TH1D* nUrqmdElMc = (TH1D*) fh_mean_nof_urqmdElectrons[kMc]->Clone();
+ TH1D* nUrqmdPosMc = (TH1D*) fh_mean_nof_urqmdPositrons[kMc]->Clone();
+ TH1D* nUrqmdElAcc = (TH1D*) fh_mean_nof_urqmdElectrons[kAcc]->Clone();
+ TH1D* nUrqmdPosAcc = (TH1D*) fh_mean_nof_urqmdPositrons[kAcc]->Clone();
+ TH1D* nUrqmdElReco = (TH1D*) fh_mean_nof_urqmdElectrons[kReco]->Clone();
+ TH1D* nUrqmdPosReco = (TH1D*) fh_mean_nof_urqmdPositrons[kReco]->Clone();
+ TH1D* nUrqmdElElid = (TH1D*) fh_mean_nof_urqmdElectrons[kElId]->Clone();
+ TH1D* nUrqmdPosElid = (TH1D*) fh_mean_nof_urqmdPositrons[kElId]->Clone();
+ TH1D* nUrqmdElTt = (TH1D*) fh_mean_nof_urqmdElectrons[kTtCut]->Clone();
+ TH1D* nUrqmdPosTt = (TH1D*) fh_mean_nof_urqmdPositrons[kTtCut]->Clone();
+ TH1D* nUrqmdElPt = (TH1D*) fh_mean_nof_urqmdElectrons[kPtCut]->Clone();
+ TH1D* nUrqmdPosPt = (TH1D*) fh_mean_nof_urqmdPositrons[kPtCut]->Clone();
+
+ int n = 2;
+ nPlutoElMc->Rebin(n * nRebin);
+ nPlutoElMc->Scale(1. / (n * nRebin));
+ nPlutoPosMc->Rebin(n * nRebin);
+ nPlutoPosMc->Scale(1. / (n * nRebin));
+ nPlutoElAcc->Rebin(n * nRebin);
+ nPlutoElAcc->Scale(1. / (n * nRebin));
+ nPlutoPosAcc->Rebin(n * nRebin);
+ nPlutoPosAcc->Scale(1. / (n * nRebin));
+ nPlutoElReco->Rebin(n * nRebin);
+ nPlutoElReco->Scale(1. / (n * nRebin));
+ nPlutoPosReco->Rebin(n * nRebin);
+ nPlutoPosReco->Scale(1. / (n * nRebin));
+ nPlutoElElid->Rebin(n * nRebin);
+ nPlutoElElid->Scale(1. / (n * nRebin));
+ nPlutoPosElid->Rebin(n * nRebin);
+ nPlutoPosElid->Scale(1. / (n * nRebin));
+ nPlutoElTt->Rebin(n * nRebin);
+ nPlutoElTt->Scale(1. / (n * nRebin));
+ nPlutoPosTt->Rebin(n * nRebin);
+ nPlutoPosTt->Scale(1. / (n * nRebin));
+ nPlutoElPt->Rebin(n * nRebin);
+ nPlutoElPt->Scale(1. / (n * nRebin));
+ nPlutoPosPt->Rebin(n * nRebin);
+ nPlutoPosPt->Scale(1. / (n * nRebin));
+
+ nUrqmdElMc->Rebin(n * nRebin);
+ nUrqmdElMc->Scale(1. / (n * nRebin));
+ nUrqmdPosMc->Rebin(n * nRebin);
+ nUrqmdPosMc->Scale(1. / (n * nRebin));
+ nUrqmdElAcc->Rebin(n * nRebin);
+ nUrqmdElAcc->Scale(1. / (n * nRebin));
+ nUrqmdPosAcc->Rebin(n * nRebin);
+ nUrqmdPosAcc->Scale(1. / (n * nRebin));
+ nUrqmdElReco->Rebin(n * nRebin);
+ nUrqmdElReco->Scale(1. / (n * nRebin));
+ nUrqmdPosReco->Rebin(n * nRebin);
+ nUrqmdPosReco->Scale(1. / (n * nRebin));
+ nUrqmdElElid->Rebin(n * nRebin);
+ nUrqmdElElid->Scale(1. / (n * nRebin));
+ nUrqmdPosElid->Rebin(n * nRebin);
+ nUrqmdPosElid->Scale(1. / (n * nRebin));
+ nUrqmdElTt->Rebin(n * nRebin);
+ nUrqmdElTt->Scale(1. / (n * nRebin));
+ nUrqmdPosTt->Rebin(n * nRebin);
+ nUrqmdPosTt->Scale(1. / (n * nRebin));
+ nUrqmdElPt->Rebin(n * nRebin);
+ nUrqmdElPt->Scale(1. / (n * nRebin));
+ nUrqmdPosPt->Rebin(n * nRebin);
+ nUrqmdPosPt->Scale(1. / (n * nRebin));
+
+ double min1 = 1e-9;
+ double max1 = 10;
+ nPlutoElMc->SetMinimum(min1);
+ nPlutoElMc->SetMaximum(max1);
+ nPlutoPosMc->SetMinimum(min1);
+ nPlutoPosMc->SetMaximum(max1);
+ nPlutoElAcc->SetMinimum(min1);
+ nPlutoElAcc->SetMaximum(max1);
+ nPlutoPosAcc->SetMinimum(min1);
+ nPlutoPosAcc->SetMaximum(max1);
+ nPlutoElReco->SetMinimum(min1);
+ nPlutoElReco->SetMaximum(max1);
+ nPlutoPosReco->SetMinimum(min1);
+ nPlutoPosReco->SetMaximum(max1);
+ nPlutoElElid->SetMinimum(min1);
+ nPlutoElElid->SetMaximum(max1);
+ nPlutoPosElid->SetMinimum(min1);
+ nPlutoPosElid->SetMaximum(max1);
+ nPlutoElTt->SetMinimum(min1);
+ nPlutoElTt->SetMaximum(max1);
+ nPlutoPosTt->SetMinimum(min1);
+ nPlutoPosTt->SetMaximum(max1);
+ nPlutoElPt->SetMinimum(min1);
+ nPlutoElPt->SetMaximum(max1);
+ nPlutoPosPt->SetMinimum(min1);
+ nPlutoPosPt->SetMaximum(max1);
+
+ nUrqmdElMc->SetMinimum(min1);
+ nUrqmdElMc->SetMaximum(max1);
+ nUrqmdPosMc->SetMinimum(min1);
+ nUrqmdPosMc->SetMaximum(max1);
+ nUrqmdElAcc->SetMinimum(min1);
+ nUrqmdElAcc->SetMaximum(max1);
+ nUrqmdPosAcc->SetMinimum(min1);
+ nUrqmdPosAcc->SetMaximum(max1);
+ nUrqmdElReco->SetMinimum(min1);
+ nUrqmdElReco->SetMaximum(max1);
+ nUrqmdPosReco->SetMinimum(min1);
+ nUrqmdPosReco->SetMaximum(max1);
+ nUrqmdElElid->SetMinimum(min1);
+ nUrqmdElElid->SetMaximum(max1);
+ nUrqmdPosElid->SetMinimum(min1);
+ nUrqmdPosElid->SetMaximum(max1);
+ nUrqmdElTt->SetMinimum(min1);
+ nUrqmdElTt->SetMaximum(max1);
+ nUrqmdPosTt->SetMinimum(min1);
+ nUrqmdPosTt->SetMaximum(max1);
+ nUrqmdElPt->SetMinimum(min1);
+ nUrqmdElPt->SetMaximum(max1);
+ nUrqmdPosPt->SetMinimum(min1);
+ nUrqmdPosPt->SetMaximum(max1);
+
+ nPlutoElMc->GetYaxis()->SetTitle("per event");
+ nPlutoElMc->SetTitle("PLUTO electrons");
+ nUrqmdElMc->GetYaxis()->SetTitle("per event");
+ nUrqmdElMc->SetTitle("UrQMD electrons");
+
+ fHM[0]->CreateCanvas("fh_nof_plutoElPos", "fh_nof_plutoElPos", 800, 800);
+ DrawH1({nPlutoElMc, nPlutoPosMc, nPlutoElAcc, nPlutoPosAcc, nPlutoElReco, nPlutoPosReco, nPlutoElElid, nPlutoPosElid,
+ nPlutoElTt, nPlutoPosTt, nPlutoElPt, nPlutoPosPt},
+ {"", "", "", "", "", "", "", "", "", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofPP = new TLegend(0.65, 0.6, 0.88, 0.93);
+ legendNofPP->SetFillColor(kWhite);
+ legendNofPP->AddEntry(nPlutoElMc, "electrons kMc");
+ legendNofPP->AddEntry(nPlutoPosMc, "positrons kMc");
+ legendNofPP->AddEntry(nPlutoElAcc, "electrons kAcc");
+ legendNofPP->AddEntry(nPlutoPosAcc, "positrons kAcc");
+ legendNofPP->AddEntry(nPlutoElReco, "electrons kReco");
+ legendNofPP->AddEntry(nPlutoPosReco, "positrons kReco");
+ legendNofPP->AddEntry(nPlutoElElid, "electrons kElId");
+ legendNofPP->AddEntry(nPlutoPosElid, "positrons kElId");
+ legendNofPP->AddEntry(nPlutoElTt, "electrons kTt");
+ legendNofPP->AddEntry(nPlutoPosTt, "positrons kTt");
+ legendNofPP->AddEntry(nPlutoElPt, "electrons kPt");
+ legendNofPP->AddEntry(nPlutoPosPt, "positrons kPt");
+ legendNofPP->Draw();
+
+ fHM[0]->CreateCanvas("fh_nof_urqmdElPos", "fh_nof_urqmdElPos", 800, 800);
+ DrawH1({nUrqmdElMc, nUrqmdPosMc, nUrqmdElAcc, nUrqmdPosAcc, nUrqmdElReco, nUrqmdPosReco, nUrqmdElElid, nUrqmdPosElid,
+ nUrqmdElTt, nUrqmdPosTt, nUrqmdElPt, nUrqmdPosPt},
+ {"", "", "", "", "", "", "", "", "", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofUP = new TLegend(0.65, 0.6, 0.88, 0.93);
+ legendNofUP->SetFillColor(kWhite);
+ legendNofUP->AddEntry(nUrqmdElMc, "electrons kMc");
+ legendNofUP->AddEntry(nUrqmdPosMc, "positrons kMc");
+ legendNofUP->AddEntry(nUrqmdElAcc, "electrons kAcc");
+ legendNofUP->AddEntry(nUrqmdPosAcc, "positrons kAcc");
+ legendNofUP->AddEntry(nUrqmdElReco, "electrons kReco");
+ legendNofUP->AddEntry(nUrqmdPosReco, "positrons kReco");
+ legendNofUP->AddEntry(nUrqmdElElid, "electrons kElId");
+ legendNofUP->AddEntry(nUrqmdPosElid, "positrons kElId");
+ legendNofUP->AddEntry(nUrqmdElTt, "electrons kTt");
+ legendNofUP->AddEntry(nUrqmdPosTt, "positrons kTt");
+ legendNofUP->AddEntry(nUrqmdElPt, "electrons kPt");
+ legendNofUP->AddEntry(nUrqmdPosPt, "positrons kPt");
+ legendNofUP->Draw();
+
+ /* 2) Draw Pair Yields */
+
+ // draw ratio e-e-/e+e+
+ TH1D* ratPairsElid = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[kElId]->Clone();
+ TH1D* ratPairsPt = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[kPtCut]->Clone();
+
+ ratPairsElid->Divide(fh_mean_combPairsPP_sameEvent_minv[kElId]);
+ ratPairsPt->Divide(fh_mean_combPairsPP_sameEvent_minv[kPtCut]);
+
+ ratPairsElid->GetXaxis()->SetRangeUser(0, 2.);
+ ratPairsElid->GetXaxis()->SetTitle(xTitle.c_str());
+ ratPairsElid->GetYaxis()->SetTitle("Ratio");
+ ratPairsElid->SetTitle("Ratio e^{-}e^{-}/e^{+}e^{+}");
+
+ fHM[0]->CreateCanvas("minv_CB_1_Ratio_eMeP", "minv_CB_1_Ratio_eMeP", 800, 800);
+ DrawH1({ratPairsElid, ratPairsPt}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "");
+
+ TLegend* legendRatio = new TLegend(0.2, 0.2, 0.7, 0.35);
+ legendRatio->SetFillColor(kWhite);
+ legendRatio->AddEntry(ratPairsElid, "ratio e^{-}e^{-}/e^{+}e^{+} pairs after El ID");
+ legendRatio->AddEntry(ratPairsPt, "ratio e^{-}e^{-}/e^{+}e^{+} pairs after Pt Cut");
+ legendRatio->Draw();
+
+ // calculate ratio of e-e-/e+e+ pairs
+ double nElPairsElid = fh_mean_combPairsMM_sameEvent_minv[kElId]->GetEntries();
+ double nPosPairsElid = fh_mean_combPairsPP_sameEvent_minv[kElId]->GetEntries();
+ double nElPairsPt = fh_mean_combPairsMM_sameEvent_minv[kPtCut]->GetEntries();
+ double nPosPairsPt = fh_mean_combPairsPP_sameEvent_minv[kPtCut]->GetEntries();
+
+ double ratioEmEpElid = nElPairsElid / nPosPairsElid;
+ double ratioEmEpPt = nElPairsPt / nPosPairsPt;
+ cout << "Ratio e-e-/e+e+ (pairs, El ID) = " << ratioEmEpElid << endl;
+ cout << "Ratio e-e-/e+e+ (pairs, Pt cut) = " << ratioEmEpPt << endl;
+
+ // draw raw pairs
+ TH1D* h21PMElidRaw = (TH1D*) fh_mean_combPairsPM_sameEvent_minv_raw[kElId]->Clone();
+ TH1D* h21PPElidRaw = (TH1D*) fh_mean_combPairsPP_sameEvent_minv_raw[kElId]->Clone();
+ TH1D* h21MMElidRaw = (TH1D*) fh_mean_combPairsMM_sameEvent_minv_raw[kElId]->Clone();
+ TH1D* h21pmElidRaw = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv_raw[kElId]->Clone();
+ TH1D* h21ppElidRaw = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv_raw[kElId]->Clone();
+ TH1D* h21mmElidRaw = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv_raw[kElId]->Clone();
+
+ h21PMElidRaw->GetXaxis()->SetTitle(xTitle.c_str());
+ h21PMElidRaw->GetYaxis()->SetTitle("absolute number");
+ h21PMElidRaw->SetTitle("same events");
+ h21pmElidRaw->GetXaxis()->SetTitle(xTitle.c_str());
+ h21pmElidRaw->GetYaxis()->SetTitle("absolute number");
+ h21pmElidRaw->SetTitle("mixed events");
+
+ fHM[0]->CreateCanvas("minv_CB_2_rawPairs_same", "minv_CB_2_rawPairs_same", 800, 800);
+ DrawH1({h21PMElidRaw, h21PPElidRaw, h21MMElidRaw}, {"", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "pe1");
+
+ TLegend* legend21 = new TLegend(0.5, 0.8, 0.92, 0.92);
+ legend21->SetFillColor(kWhite);
+ legend21->AddEntry(h21PMElidRaw, "e+e- pairs after El ID");
+ legend21->AddEntry(h21PPElidRaw, "e+e+ pairs after El ID");
+ legend21->AddEntry(h21MMElidRaw, "e-e- pairs after El ID");
+ legend21->Draw();
+
+ fHM[0]->CreateCanvas("minv_CB_2_rawPairs_mixed", "minv_CB_2_rawPairs_mixed", 800, 800);
+ DrawH1({h21pmElidRaw, h21ppElidRaw, h21mmElidRaw}, {"", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "pe1");
+
+ TLegend* legend21b = new TLegend(0.5, 0.8, 0.92, 0.92);
+ legend21b->SetFillColor(kWhite);
+ legend21b->AddEntry(h21pmElidRaw, "e+e- pairs after El ID");
+ legend21b->AddEntry(h21ppElidRaw, "e+e+ pairs after El ID");
+ legend21b->AddEntry(h21mmElidRaw, "e-e- pairs after El ID");
+ legend21b->Draw();
+
+ // draw scaled pairs (same event)
+ TH1D* Bpm = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone();
+ TH1D* Bpp = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[kElId]->Clone();
+ TH1D* Bmm = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[kElId]->Clone();
+
+ double bW = Bpm->GetBinWidth(1); // MIND: 'bW' is used throughout this draw method!
+ cout << "DrawMinvCombSignalAndBg(): bW = " << bW << endl;
+
+ Bpm->GetXaxis()->SetRangeUser(0, 2.);
+ Bpm->GetXaxis()->SetTitle(xTitle.c_str());
+ Bpm->GetYaxis()->SetTitle(yTitle.c_str());
+ Bpm->SetTitle("same event");
+
+ fHM[0]->CreateCanvas("minv_CB_2_pairsSameEvent", "minv_CB_2_pairsSameEvent", 800, 800);
+ DrawH1({Bpm, Bpp, Bmm}, {"", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendB = new TLegend(0.5, 0.75, 0.95, 0.93);
+ legendB->SetFillColor(kWhite);
+ legendB->AddEntry(Bpm, "e+e- pairs (Elid)");
+ legendB->AddEntry(Bpp, "e+e+ pairs (Elid)");
+ legendB->AddEntry(Bmm, "e-e- pairs (Elid)");
+ legendB->Draw();
+
+ // draw scaled pairs (mixed events)
+ TH1D* bpm = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv[kElId]->Clone();
+ TH1D* bpp = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv[kElId]->Clone();
+ TH1D* bmm = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv[kElId]->Clone();
+
+ bpm->GetXaxis()->SetRangeUser(0, 2.);
+ bpm->GetXaxis()->SetTitle(xTitle.c_str());
+ bpm->GetYaxis()->SetTitle(yTitle.c_str());
+ bpm->SetTitle("mixed events");
+
+ fHM[0]->CreateCanvas("minv_CB_2_pairsMixedEvents", "minv_CB_2_pairsMixedEvents", 800, 800);
+ DrawH1({bpm, bpp, bmm}, {"", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendb = new TLegend(0.5, 0.75, 0.95, 0.93);
+ legendb->SetFillColor(kWhite);
+ legendb->AddEntry(bpm, "e+e- pairs (Elid)");
+ legendb->AddEntry(bpp, "e+e+ pairs (Elid)");
+ legendb->AddEntry(bmm, "e-e- pairs (Elid)");
+ legendb->Draw();
+
+ // compare B++ with b++ (and -- and +-, resp.); therefor normalize to integral(400,700)
+ // convert MeV into Bin
+ int intFrom = 400; // lower and upper range to normalize (in Mev/c^2)
+ int intTo = 700;
+ int nBinsOrig = 4000; // was set in CbmAnaDielectronTask.cxx
+ int upperLimit = 4000; // in MeV (assumed, that lower limit = 0)
+ double bin400 = (double) (1. * nBinsOrig / nRebin) * (1. * intFrom / upperLimit);
+ double bin700 = (double) (1. * nBinsOrig / nRebin) * (1. * intTo / upperLimit);
+ cout << "bin400 = " << bin400 << endl;
+ cout << "bin700 = " << bin700 << endl;
+
+ // e+e- pairs
+ TH1D* BpmNormed = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone(); // norm. factor = 1
+ TH1D* bpmNormed = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv[kElId]->Clone();
+
+ BpmNormed->GetXaxis()->SetRangeUser(0, 2.);
+ BpmNormed->GetXaxis()->SetTitle(xTitle.c_str());
+ BpmNormed->GetYaxis()->SetTitle(yTitle.c_str());
+
+ double intSamePM = BpmNormed->Integral(bin400, bin700);
+ double intMixedPM = bpmNormed->Integral(bin400, bin700);
+ double normPM = (double) intSamePM / intMixedPM; // mixed histograms will be scaled with this factor
+
+ int nofBins = BpmNormed->GetNbinsX(); // MIND: 'nofBins' is used throughout this draw method!
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double content = bpmNormed->GetBinContent(iBin);
+ content *= normPM;
+ bpmNormed->SetBinContent(iBin, content);
+ }
- gPad->SetLogy(true);
+ // e+e+ pairs
+ TH1D* BppNormed = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[kElId]->Clone(); // norm. factor = 1
+ TH1D* bppNormed = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv[kElId]->Clone();
+
+ BppNormed->GetXaxis()->SetRangeUser(0, 2.);
+ BppNormed->GetXaxis()->SetTitle(xTitle.c_str());
+ BppNormed->GetYaxis()->SetTitle(yTitle.c_str());
+
+ double intSamePP = BppNormed->Integral(bin400, bin700);
+ double intMixedPP = bppNormed->Integral(bin400, bin700);
+ double normPP = (double) intSamePP / intMixedPP; // mixed histograms will be scaled with this factor
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double content = bppNormed->GetBinContent(iBin);
+ content *= normPP;
+ bppNormed->SetBinContent(iBin, content);
+ }
+
+ // e-e- pairs
+ TH1D* BmmNormed = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[kElId]->Clone(); // norm. factor = 1
+ TH1D* bmmNormed = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv[kElId]->Clone();
+
+ BmmNormed->GetXaxis()->SetRangeUser(0, 2.);
+ BmmNormed->GetXaxis()->SetTitle(xTitle.c_str());
+ BmmNormed->GetYaxis()->SetTitle(yTitle.c_str());
+
+ double intSameMM = BmmNormed->Integral(bin400, bin700);
+ double intMixedMM = bmmNormed->Integral(bin400, bin700);
+ double normMM = (double) intSameMM / intMixedMM; // mixed histograms will be scaled with this factor
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double content = bmmNormed->GetBinContent(iBin);
+ content *= normMM;
+ bmmNormed->SetBinContent(iBin, content);
+ }
+
+ BpmNormed->SetTitle("compare same/mixed e+e- events");
+ BppNormed->SetTitle("compare same/mixed e+e+ events");
+ BmmNormed->SetTitle("compare same/mixed e-e- events");
+
+ TCanvas* c = fHM[0]->CreateCanvas("minv_CB_2_compare_same-mixed", "minv_CB_2_compare_same-mixed", 1800, 600);
+ c->Divide(3, 1);
+
+ c->cd(1);
+ if (setMinMax) {
+ BpmNormed->SetMinimum(yMin);
+ BpmNormed->SetMaximum(2 * yMax);
+ }
+ DrawH1({BpmNormed, bpmNormed}, {"p", "p"}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+ TLegend* leg1 = new TLegend(0.5, 0.78, 0.92, 0.9);
+ leg1->SetFillColor(kWhite);
+ leg1->AddEntry(BpmNormed, "e^{+}-e^{-} same event (Elid)");
+ leg1->AddEntry(bpmNormed, "e^{+}-e^{-} mixed events (Elid)");
+ leg1->Draw();
+
+ c->cd(2);
+ if (setMinMax) {
+ BppNormed->SetMinimum(yMin);
+ BppNormed->SetMaximum(2 * yMax);
+ }
+ DrawH1({BppNormed, bppNormed}, {"p", "p"}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+ TLegend* leg2 = new TLegend(0.5, 0.78, 0.92, 0.9);
+ leg2->SetFillColor(kWhite);
+ leg2->AddEntry(BppNormed, "e^{+}-e^{+} same event (Elid)");
+ leg2->AddEntry(bppNormed, "e^{+}-e^{+} mixed events (Elid)");
+ leg2->Draw();
+
+ c->cd(3);
+ if (setMinMax) {
+ BmmNormed->SetMinimum(yMin);
+ BmmNormed->SetMaximum(2 * yMax);
+ }
+ DrawH1({BmmNormed, bmmNormed}, {"p", "p"}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+ TLegend* leg3 = new TLegend(0.5, 0.78, 0.92, 0.9);
+ leg3->SetFillColor(kWhite);
+ leg3->AddEntry(BmmNormed, "e^{-}-e^{-} same event (Elid)");
+ leg3->AddEntry(bmmNormed, "e^{-}-e^{-} mixed events (Elid)");
+ leg3->Draw();
+
+ // compare B+- with (Cocktail + BG)
+ TH1D* Bpm1Elid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId];
+ TH1D* Bpm1Pt = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kPtCut];
+ TH1D* cocBg1Elid = (TH1D*) GetCoctailMinv(kElId);
+ TH1D* cocBg1Pt = (TH1D*) GetCoctailMinv(kPtCut);
+
+ cocBg1Elid->Add(fh_mean_bg_minv[kElId]);
+ cocBg1Pt->Add(fh_mean_bg_minv[kPtCut]);
+
+ cocBg1Elid->Rebin(nRebin);
+ cocBg1Elid->Scale(1. / bW);
+ cocBg1Pt->Rebin(nRebin);
+ cocBg1Pt->Scale(1. / bW);
+
+ Bpm1Elid->GetXaxis()->SetRangeUser(0, 2.);
+ Bpm1Elid->GetXaxis()->SetTitle(xTitle.c_str());
+ Bpm1Elid->GetYaxis()->SetTitle(yTitle.c_str());
+
+ fHM[0]->CreateCanvas("minv_CB_2_BpmVsMC", "minv_CB_2_BpmVsMC", 800, 800);
+ DrawH1({Bpm1Elid, cocBg1Elid, Bpm1Pt, cocBg1Pt}, {"", "", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99,
+ "HIST L");
+
+ TLegend* legendCocBg = new TLegend(0.5, 0.7, 0.92, 0.91);
+ legendCocBg->SetFillColor(kWhite);
+ legendCocBg->AddEntry(Bpm1Elid, "N_{same}^{+-} (Elid)");
+ legendCocBg->AddEntry(cocBg1Elid, "Cocktail + BG (Elid)");
+ legendCocBg->AddEntry(Bpm1Pt, "N_{same}^{+-} (P_{t} cut)");
+ legendCocBg->AddEntry(cocBg1Pt, "Cocktail + BG (P_{t} cut)");
+ legendCocBg->Draw();
+
+ TH1D* ratCocBgBpmElid = (TH1D*) cocBg1Elid->Clone();
+ TH1D* ratCocBgBpmPt = (TH1D*) cocBg1Pt->Clone();
+
+ ratCocBgBpmElid->Divide(Bpm1Elid);
+ ratCocBgBpmPt->Divide(Bpm1Pt);
+
+ ratCocBgBpmElid->GetXaxis()->SetRangeUser(0, 2.);
+ ratCocBgBpmElid->GetYaxis()->SetTitle(0);
+ ratCocBgBpmElid->SetTitle("Ratio (Cocktail + BG) / N_{same}^{+-}");
+
+ fHM[0]->CreateCanvas("minv_CB_2_BpmVsMC2", "minv_CB_2_BpmVsMC2", 800, 800);
+ DrawH1({ratCocBgBpmElid, ratCocBgBpmPt}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendCocBg2 = new TLegend(0.55, 0.75, 0.95, 0.91);
+ legendCocBg2->SetFillColor(kWhite);
+ legendCocBg2->AddEntry(ratCocBgBpmElid, "El ID cut");
+ legendCocBg2->AddEntry(ratCocBgBpmPt, "P_{t} cut");
+ legendCocBg2->Draw();
+
+ /* 3) Draw Geometric Mean */
+ TH1D* geomMeanSameElid = (TH1D*) fh_mean_combBg_GeomMeanSame_minv[kElId]->Clone();
+ TH1D* geomMeanMixedElid = (TH1D*) fh_mean_combBg_GeomMeanMixed_minv[kElId]->Clone();
+ TH1D* geomMeanSamePt = (TH1D*) fh_mean_combBg_GeomMeanSame_minv[kPtCut]->Clone();
+ TH1D* geomMeanMixedPt = (TH1D*) fh_mean_combBg_GeomMeanMixed_minv[kPtCut]->Clone();
+
+ double intGMsameElid = geomMeanSameElid->Integral(bin400, bin700);
+ double intGMmixedElid = geomMeanMixedElid->Integral(bin400, bin700);
+ double normGMElid = (double) intGMsameElid / intGMmixedElid; // mixed histograms will be scaled with this factor
+
+ double intGMsamePt = geomMeanSamePt->Integral(bin400, bin700);
+ double intGMmixedPt = geomMeanMixedPt->Integral(bin400, bin700);
+ double normGMPt = (double) intGMsamePt / intGMmixedPt;
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double contentElid = geomMeanMixedElid->GetBinContent(iBin);
+ double contentPt = geomMeanMixedPt->GetBinContent(iBin);
+ contentElid *= normGMElid;
+ contentPt *= normGMPt;
+ geomMeanMixedElid->SetBinContent(iBin, contentElid);
+ geomMeanMixedPt->SetBinContent(iBin, contentPt);
+ }
+ if (setMinMax) {
+ geomMeanSameElid->SetMinimum(yMin);
+ geomMeanSameElid->SetMaximum(yMax);
+ }
+
+ geomMeanSameElid->GetXaxis()->SetRangeUser(0, 2.);
+ geomMeanSameElid->GetXaxis()->SetTitle(xTitle.c_str());
+ geomMeanSameElid->GetYaxis()->SetTitle(yTitle.c_str());
+ geomMeanSameElid->SetTitle("geometric mean (normalized)");
+ geomMeanMixedElid->SetTitle(0);
+
+ fHM[0]->CreateCanvas("minv_CB_3_geomMean", "minv_CB_3_geomMean", 800, 800);
+ DrawH1({geomMeanSameElid, geomMeanMixedElid, geomMeanSamePt, geomMeanMixedPt}, {"", "", "", ""}, kLinear, kLog, false,
+ 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendgeomMean = new TLegend(0.55, 0.8, 0.92, 0.91);
+ legendgeomMean->SetFillColor(kWhite);
+ legendgeomMean->AddEntry(geomMeanSameElid, "same event (elid)");
+ legendgeomMean->AddEntry(geomMeanMixedElid, "mixed events (elid)");
+ legendgeomMean->AddEntry(geomMeanSamePt, "same event (Pt)");
+ legendgeomMean->AddEntry(geomMeanMixedPt, "mixed events (Pt)");
+ legendgeomMean->Draw();
+
+ //draw geom. mean, assembled of same (<= 0.3 GeV) and normed mixed (> 0.3 GeV) data
+
+ // convert MeV into Bin
+ int rFrom = 300; // lower and upper range to normalize (in Mev/c^2)
+ int rTo = 1000;
+ double binFrom = (double) (1. * nBinsOrig / nRebin) * (1. * rFrom / upperLimit);
+ double binTo = (double) (1. * nBinsOrig / nRebin) * (1. * rTo / upperLimit);
+
+ TH1D* geomMeanAss = (TH1D*) fh_mean_combBg_GeomMeanSame_minv[kElId]->Clone();
+
+ double intSameGMsp = geomMeanAss->Integral(binFrom, binTo);
+ double intMixedGMsp = geomMeanMixedElid->Integral(binFrom, binTo);
+ double normGMsp = (double) intSameGMsp / intMixedGMsp;
+
+ for (int iBin = binFrom + 1; iBin <= nofBins; iBin++) { // from 300 MeV on normalized data from mixed event
+ double cont = geomMeanMixedElid->GetBinContent(iBin);
+ cont *= normGMsp;
+ geomMeanAss->SetBinContent(iBin, cont);
+ }
+ if (setMinMax) {
+ geomMeanAss->SetMinimum(yMin);
+ geomMeanAss->SetMaximum(yMax);
+ }
+
+ geomMeanAss->GetXaxis()->SetRangeUser(0, 2.);
+ geomMeanAss->GetXaxis()->SetTitle(xTitle.c_str());
+ geomMeanAss->GetYaxis()->SetTitle(yTitle.c_str());
+ geomMeanAss->SetTitle("geometric mean (splitted)");
+
+ fHM[0]->CreateCanvas("minv_CB_3_geomMean_assembled", "minv_CB_3_geomMean_assembled", 800, 800);
+ DrawH1({geomMeanAss}, {""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist l");
+
+ TLegend* legendgeomMeanAss = new TLegend(0.55, 0.8, 0.91, 0.91);
+ legendgeomMeanAss->SetFillColor(kWhite);
+ legendgeomMeanAss->AddEntry(geomMeanAss, "geom. mean");
+ legendgeomMeanAss->Draw();
+
+ /* 4) Draw B_comb with Coc+BG and with Bpm */
+ TH1D* cbgElid = (TH1D*) fh_mean_combBg_minv[kElId]->Clone();
+ TH1D* cbgPtCut = (TH1D*) fh_mean_combBg_minv[kPtCut]->Clone();
+
+ TH1D* bgCocElid = (TH1D*) GetCoctailMinv(kElId);
+ TH1D* bgCocPt = (TH1D*) GetCoctailMinv(kPtCut);
+
+ bgCocElid->Add(fh_mean_bg_minv[kElId]);
+ bgCocPt->Add(fh_mean_bg_minv[kPtCut]);
+
+ TH1D* BpmElid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone(); // MIND: BpmElid is also used below!
+ TH1D* BpmPt = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kPtCut]->Clone();
+
+ bgCocElid->Rebin(nRebin); // these histos have to be rebinned and scaled here because not done yet as with CB histos
+ bgCocPt->Rebin(nRebin);
+ bgCocElid->Scale(1. / bW);
+ bgCocPt->Scale(1. / bW);
+
+ cbgElid->GetXaxis()->SetRangeUser(0, 2.);
+ cbgElid->GetXaxis()->SetTitle(xTitle.c_str());
+ cbgElid->GetYaxis()->SetTitle(yTitle.c_str());
+ cbgElid->SetTitle(0);
+
+ if (setMinMax) {
+ cbgElid->SetMinimum(yMin);
+ cbgElid->SetMaximum(yMax);
+ }
+
+ fHM[0]->CreateCanvas("minv_CB_4_bComb", "minv_CB_4_bComb", 800, 800);
+ DrawH1({cbgElid, bgCocElid}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendBg = new TLegend(0.55, 0.7, 0.95, 0.9);
+ legendBg->SetFillColor(kWhite);
+ legendBg->AddEntry(cbgElid, "comb. BG (2 * k * geomMean, elid)");
+ legendBg->AddEntry(bgCocElid, "Cocktail + BG (elid)");
+ legendBg->Draw();
+
+ TH1D* ratCocBgCbgElid = (TH1D*) bgCocElid->Clone();
+ TH1D* ratCocBgCbgPt = (TH1D*) bgCocPt->Clone();
+
+ ratCocBgCbgElid->Divide(cbgElid);
+ ratCocBgCbgPt->Divide(cbgPtCut);
+
+ if (setMinMax) {
+ ratCocBgCbgElid->SetMinimum(2e-1);
+ ratCocBgCbgElid->SetMaximum(50);
+ ratCocBgCbgPt->SetMinimum(2e-1);
+ ratCocBgCbgPt->SetMaximum(50);
+ cbgElid->SetMinimum(yMin);
+ cbgElid->SetMaximum(yMax);
+ }
+
+ fHM[0]->CreateCanvas("minv_CB_4_bCombVsBpm", "minv_CB_4_bCombVsBpm", 800, 800);
+ DrawH1({cbgElid, BpmElid, cbgPtCut, BpmPt}, {"", "", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendBg2 = new TLegend(0.55, 0.7, 0.92, 0.91);
+ legendBg2->SetFillColor(kWhite);
+ legendBg2->AddEntry(cbgElid, "comb. BG, Elid");
+ legendBg2->AddEntry(BpmElid, "N_{same}^{+-}, Elid");
+ legendBg2->AddEntry(cbgPtCut, "comb. BG, PtCut");
+ legendBg2->AddEntry(BpmPt, "N_{same}^{+-}, PtCut");
+ legendBg2->Draw();
+
+ ratCocBgCbgElid->GetXaxis()->SetRangeUser(0, 2.);
+ ratCocBgCbgElid->SetTitle("Ratio (Cocktail + BG) / B_{c}");
+
+ fHM[0]->CreateCanvas("minv_CB_4_bComb_rat", "minv_CB_4_bComb_rat", 800, 800);
+ DrawH1({ratCocBgCbgElid, ratCocBgCbgPt}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendBg1 = new TLegend(0.55, 0.8, 0.92, 0.9);
+ legendBg1->SetFillColor(kWhite);
+ legendBg1->AddEntry(ratCocBgCbgElid, "El ID cut");
+ legendBg1->AddEntry(ratCocBgCbgPt, "Pt cut");
+ legendBg1->Draw();
+
+ // draw ass. CB with all steps
+ TH1D* cbgAssElid = (TH1D*) fh_mean_combBg_assemb_minv[kElId]->Clone();
+ TH1D* cbgAssGamma = (TH1D*) fh_mean_combBg_assemb_minv[kGammaCut]->Clone();
+ TH1D* cbgAssMvd1 = (TH1D*) fh_mean_combBg_assemb_minv[kMvd1Cut]->Clone();
+ TH1D* cbgAssSt = (TH1D*) fh_mean_combBg_assemb_minv[kStCut]->Clone();
+ TH1D* cbgAssRt = (TH1D*) fh_mean_combBg_assemb_minv[kRtCut]->Clone();
+ TH1D* cbgAssTt = (TH1D*) fh_mean_combBg_assemb_minv[kTtCut]->Clone();
+ TH1D* cbgAssPt = (TH1D*) fh_mean_combBg_assemb_minv[kPtCut]->Clone();
+
+ TH1D* ratCocBgCbgAssElid = (TH1D*) bgCocElid->Clone();
+ TH1D* ratCocBgCbgAssPt = (TH1D*) bgCocPt->Clone();
+
+ ratCocBgCbgAssElid->Divide(cbgAssElid);
+ ratCocBgCbgAssPt->Divide(cbgAssPt);
+
+ if (setMinMax) {
+ ratCocBgCbgAssElid->SetMinimum(2e-1);
+ ratCocBgCbgAssElid->SetMaximum(50);
+ ratCocBgCbgAssPt->SetMinimum(2e-1);
+ ratCocBgCbgAssPt->SetMaximum(50);
+ }
+
+ ratCocBgCbgAssElid->GetXaxis()->SetRangeUser(0, 2.);
+ ratCocBgCbgAssElid->SetTitle("Ratio (Cocktail + BG) / B_{c_ass}");
+
+ fHM[0]->CreateCanvas("minv_CB_4_bComb_rat2", "minv_CB_4_bComb_rat2", 800, 800);
+ DrawH1({ratCocBgCbgAssElid, ratCocBgCbgAssPt}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendBg4 = new TLegend(0.55, 0.8, 0.92, 0.9);
+ legendBg4->SetFillColor(kWhite);
+ legendBg4->AddEntry(ratCocBgCbgAssElid, "El ID cut");
+ legendBg4->AddEntry(ratCocBgCbgAssPt, "Pt cut");
+ legendBg4->Draw();
+
+ // conv. vs assembled CB
+ TH1D* h44cbgElid = (TH1D*) fh_mean_combBg_minv[kElId]->Clone();
+ TH1D* h44cbgAssElid = (TH1D*) fh_mean_combBg_assemb_minv[kElId]->Clone();
+
+ if (setMinMax) {
+ h44cbgElid->SetMinimum(yMin);
+ h44cbgElid->SetMaximum(yMax);
+ }
+
+ h44cbgElid->GetXaxis()->SetRangeUser(0, 2.);
+ h44cbgElid->GetXaxis()->SetTitle(xTitle.c_str());
+ h44cbgElid->GetYaxis()->SetTitle(yTitle.c_str());
+ h44cbgElid->SetTitle("Combinatorial BG");
+
+ fHM[0]->CreateCanvas("minv_CB_4_bComb_convVsAss", "minv_CB_4_bComb_convVsAss", 800, 800);
+ DrawH1({h44cbgElid, h44cbgAssElid}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "pe1");
+
+ TLegend* legendBg3 = new TLegend(0.55, 0.7, 0.92, 0.91);
+ legendBg3->SetFillColor(kWhite);
+ legendBg3->AddEntry(h44cbgElid, "conv. CB, Elid");
+ legendBg3->AddEntry(h44cbgAssElid, "assembled CB, Elid");
+ legendBg3->Draw();
+
+ // assembled CB with all steps
+ if (setMinMax) {
+ cbgAssElid->SetMinimum(yMin);
+ cbgAssElid->SetMaximum(yMax);
+ }
+ cbgAssElid->GetXaxis()->SetRangeUser(0, 2.);
+ cbgAssElid->GetXaxis()->SetTitle(xTitle.c_str());
+ cbgAssElid->GetYaxis()->SetTitle(yTitle.c_str());
+ cbgAssElid->SetTitle("Combinatorial BG");
+
+ fHM[0]->CreateCanvas("minv_CB_4_bComb_assembled_steps", "minv_CB_4_bComb_assembled_steps", 800, 800);
+ DrawH1({cbgAssElid, cbgAssGamma, cbgAssMvd1, cbgAssSt, cbgAssRt, cbgAssTt, cbgAssPt}, {"", "", "", "", "", "", ""},
+ kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "pe1");
+
+ TLegend* legendBg5 = new TLegend(0.7, 0.65, 0.92, 0.91);
+ legendBg5->SetFillColor(kWhite);
+ legendBg5->AddEntry(cbgAssElid, "elid");
+ legendBg5->AddEntry(cbgAssGamma, "gamma cut");
+ legendBg5->AddEntry(cbgAssMvd1, "MVD1 cut");
+ legendBg5->AddEntry(cbgAssSt, "ST cut");
+ legendBg5->AddEntry(cbgAssRt, "RT cut");
+ legendBg5->AddEntry(cbgAssTt, "TT cut");
+ legendBg5->AddEntry(cbgAssPt, "P_{t} cut");
+ legendBg5->Draw();
+
+ /* 5) Draw k Factor */
+ TH1D* kReconst = (TH1D*) fh_mean_combBg_k_minv[kReco]->Clone();
+ TH1D* kChi2 = (TH1D*) fh_mean_combBg_k_minv[kChi2Prim]->Clone();
+ TH1D* kElid = (TH1D*) fh_mean_combBg_k_minv[kElId]->Clone();
+ TH1D* kGamma = (TH1D*) fh_mean_combBg_k_minv[kGammaCut]->Clone();
+ TH1D* kMvd1 = (TH1D*) fh_mean_combBg_k_minv[kMvd1Cut]->Clone();
+ TH1D* kSt = (TH1D*) fh_mean_combBg_k_minv[kStCut]->Clone();
+ TH1D* kRt = (TH1D*) fh_mean_combBg_k_minv[kRtCut]->Clone();
+ TH1D* kTt = (TH1D*) fh_mean_combBg_k_minv[kTtCut]->Clone();
+ TH1D* kPt = (TH1D*) fh_mean_combBg_k_minv[kPtCut]->Clone();
+
+ if (setMinMax) {
+ kReconst->SetMinimum(0.7);
+ kReconst->SetMaximum(1.3);
+ kChi2->SetMinimum(0.7);
+ kChi2->SetMaximum(1.3);
+ kElid->SetMinimum(0.7);
+ kElid->SetMaximum(1.3);
+ kGamma->SetMinimum(0.7);
+ kGamma->SetMaximum(1.3);
+ kMvd1->SetMinimum(0.7);
+ kMvd1->SetMaximum(1.3);
+ kSt->SetMinimum(0.7);
+ kSt->SetMaximum(1.3);
+ kRt->SetMinimum(0.7);
+ kRt->SetMaximum(1.3);
+ kTt->SetMinimum(0.7);
+ kTt->SetMaximum(1.3);
+ kPt->SetMinimum(0.7);
+ kPt->SetMaximum(1.3);
+ }
+
+ kReconst->GetXaxis()->SetRangeUser(0, 2.);
+ kReconst->GetXaxis()->SetTitle(xTitle.c_str());
+ kReconst->GetYaxis()->SetTitle(0);
+ kReconst->SetTitle("k factor");
+
+ fHM[0]->CreateCanvas("minv_CB_5_k", "minv_CB_5_k", 800, 800);
+ DrawH1({kReconst, kChi2, kElid, kGamma, kMvd1, kSt, kRt, kTt, kPt}, {"", "", "", "", "", "", "", "", ""}, kLinear,
+ kLinear, false, 0.8, 0.8, 0.99, 0.99, "HIST L");
+
+ TLegend* legendK = new TLegend(0.7, 0.7, 0.92, 0.91);
+ legendK->SetFillColor(kWhite);
+ legendK->AddEntry(kReconst, "reco");
+ legendK->AddEntry(kChi2, "#chi^{2} ");
+ legendK->AddEntry(kElid, "elid");
+ legendK->AddEntry(kGamma, "gamma cut");
+ legendK->AddEntry(kMvd1, "MVD1 cut");
+ legendK->AddEntry(kSt, "ST cut");
+ legendK->AddEntry(kRt, "RT cut");
+ legendK->AddEntry(kTt, "TT cut");
+ legendK->AddEntry(kPt, "P_{t} cut");
+ legendK->Draw();
+
+ /* 6) Draw Combinatorial Signal (all from assembled comb. BG) */
+ // from 'N+-same'
+
+ // elid, draw also components
+ TH1D* h61NpmElid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone();
+ TH1D* h61cbAssElid = (TH1D*) fh_mean_combBg_assemb_minv[kElId]->Clone();
+ TH1D* h61SigNpmElid = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kElId]->Clone();
+ TH1D* h61coc1Elid = (TH1D*) GetCoctailMinv(kElId);
+
+ h61coc1Elid->Add(fh_mean_bg_minv[kElId]);
+ h61coc1Elid->Add(fh_mean_bg_minv[kElId], -1.);
+
+ h61coc1Elid->Rebin(nRebin);
+ h61coc1Elid->Scale(1. / bW);
+
+ if (setMinMax) {
+ h61NpmElid->SetMinimum(yMin);
+ h61NpmElid->SetMaximum(yMax);
+ h61cbAssElid->SetMinimum(yMin);
+ h61cbAssElid->SetMaximum(yMax);
+ h61SigNpmElid->SetMinimum(yMin);
+ h61SigNpmElid->SetMaximum(yMax);
+ h61coc1Elid->SetMinimum(yMin);
+ h61coc1Elid->SetMaximum(yMax);
+ }
+
+ h61NpmElid->GetXaxis()->SetRangeUser(0, 2.);
+ h61NpmElid->GetYaxis()->SetTitle(yTitle.c_str());
+ h61NpmElid->SetTitle(0);
+ h61cbAssElid->GetXaxis()->SetRangeUser(0, 2.);
+ h61SigNpmElid->GetXaxis()->SetRangeUser(0, 2.);
+ h61coc1Elid->GetXaxis()->SetRangeUser(0, 2.);
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_Npm_elid", "minv_CB_6_signal_Npm_elid", 800, 800);
+ DrawH1({h61NpmElid, h61cbAssElid, h61SigNpmElid, h61coc1Elid}, {"", "", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99,
+ 0.99, "hist p");
+
+ TLegend* legend61 = new TLegend(0.55, 0.7, 0.98, 0.9);
+ legend61->SetFillColor(kWhite);
+ legend61->AddEntry(h61NpmElid, "N_{same}^{+-} (Elid)");
+ legend61->AddEntry(h61cbAssElid, "B_{c} (assembled, Elid)");
+ legend61->AddEntry(h61SigNpmElid, "Signal (N_{same}^{+-} - B_{c}, Elid)");
+ legend61->AddEntry(h61coc1Elid, "Cocktail, Elid)");
+ legend61->Draw();
+
+ // all steps, draw with and wo error
+ TH1D* h63SigNpmElid = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kElId]->Clone();
+ TH1D* h63SigNpmGamma = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kGammaCut]->Clone();
+ TH1D* h63SigNpmMvd1 = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kMvd1Cut]->Clone();
+ TH1D* h63SigNpmSt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kStCut]->Clone();
+ TH1D* h63SigNpmRt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kRtCut]->Clone();
+ TH1D* h63SigNpmTt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kTtCut]->Clone();
+ TH1D* h63SigNpmPt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kPtCut]->Clone();
+
+ if (setMinMax) {
+ h63SigNpmElid->SetMinimum(yMin / 10.); // setting min/max somehow makes the histos empty
+ h63SigNpmElid->SetMaximum(yMax);
+ h63SigNpmGamma->SetMinimum(yMin / 10.);
+ h63SigNpmGamma->SetMaximum(yMax);
+ h63SigNpmMvd1->SetMinimum(yMin / 10.);
+ h63SigNpmMvd1->SetMaximum(yMax);
+ h63SigNpmSt->SetMinimum(yMin / 10.);
+ h63SigNpmSt->SetMaximum(yMax);
+ h63SigNpmRt->SetMinimum(yMin / 10.);
+ h63SigNpmRt->SetMaximum(yMax);
+ h63SigNpmTt->SetMinimum(yMin / 10.);
+ h63SigNpmTt->SetMaximum(yMax);
+ h63SigNpmPt->SetMinimum(yMin / 10.);
+ h63SigNpmPt->SetMaximum(yMax);
+ }
+
+ h63SigNpmElid->GetXaxis()->SetRangeUser(0, 2.);
+ h63SigNpmElid->GetXaxis()->SetTitle(xTitle.c_str());
+ h63SigNpmElid->GetYaxis()->SetTitle(yTitle.c_str());
+ h63SigNpmElid->SetTitle("Signal");
+ /*h63SigNpmPt->GetXaxis()->SetRangeUser(0, 2.);
+ h63SigNpmPt->GetXaxis()->SetTitle(xTitle.c_str());
+ h63SigNpmPt->GetYaxis()->SetTitle(yTitle.c_str());
+ h63SigNpmPt->SetTitle("Signal");*/
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_Npm_steps", "minv_CB_6_signal_Npm_steps", 800, 800);
+ DrawH1({h63SigNpmElid, h63SigNpmGamma, h63SigNpmMvd1, h63SigNpmSt, h63SigNpmRt, h63SigNpmTt, h63SigNpmPt},
+ {"", "", "", "", "", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legend63 = new TLegend(0.7, 0.6, 0.92, 0.9);
+ legend63->SetFillColor(kWhite);
+ legend63->AddEntry(h63SigNpmElid, "elid");
+ legend63->AddEntry(h63SigNpmGamma, "gamma cut");
+ legend63->AddEntry(h63SigNpmMvd1, "MVD1 cut");
+ legend63->AddEntry(h63SigNpmSt, "ST cut");
+ legend63->AddEntry(h63SigNpmRt, "Rt cut");
+ legend63->AddEntry(h63SigNpmTt, "Tt cut");
+ legend63->AddEntry(h63SigNpmPt, "P_{t} cut");
+ legend63->Draw();
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_Npm_steps2", "minv_CB_6_signal_Npm_steps2", 800, 800);
+ DrawH1({h63SigNpmElid, h63SigNpmGamma, h63SigNpmMvd1, h63SigNpmSt, h63SigNpmRt, h63SigNpmTt, h63SigNpmPt},
+ {"", "", "", "", "", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist pe1");
+
+ TLegend* legend63e = new TLegend(0.7, 0.6, 0.92, 0.9);
+ legend63e->SetFillColor(kWhite);
+ legend63e->AddEntry(h63SigNpmElid, "elid");
+ legend63e->AddEntry(h63SigNpmGamma, "gamma cut");
+ legend63e->AddEntry(h63SigNpmMvd1, "MVD1 cut");
+ legend63e->AddEntry(h63SigNpmSt, "ST cut");
+ legend63e->AddEntry(h63SigNpmRt, "Rt cut");
+ legend63e->AddEntry(h63SigNpmTt, "Tt cut");
+ legend63e->AddEntry(h63SigNpmPt, "P_{t} cut");
+ legend63e->Draw();
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_Npm_steps3", "minv_CB_6_signal_Npm_steps3", 800, 800);
+ DrawH1({h63SigNpmElid}, {""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist pe1");
+
+ TLegend* legend63e3 = new TLegend(0.7, 0.8, 0.92, 0.9);
+ legend63e3->SetFillColor(kWhite);
+ legend63e3->AddEntry(h63SigNpmElid, "elid");
+ legend63e3->Draw();
+
+ /*fHM[0]->CreateCanvas("minv_CB_6_signal_Npm_steps4", "minv_CB_6_signal_Npm_steps4", 800, 800); Wenn aktiviert, sind Elid- und Pt-Kurven oben beide gleichfarbig
+ DrawH1({h63SigNpmPt}, {""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist pe1");
+
+ TLegend* legend63e4 = new TLegend(0.7, 0.8, 0.92, 0.9);
+ legend63e4->SetFillColor(kWhite);
+ legend63e4->AddEntry(h63SigNpmPt, "P_{t} cut");
+ legend63e4->Draw();*/
+
+ // from 'Cocktail + BG'
+ // elid
+ TH1D* h65SigCocBgElid = (TH1D*) fh_mean_combSignalBCoc_assemb_minv[kElId]->Clone();
+ TH1D* h65cbAssElid = (TH1D*) fh_mean_combBg_assemb_minv[kElId]->Clone();
+ TH1D* h65cocBgElid = (TH1D*) GetCoctailMinv(kElId);
+ TH1D* h65coc1Elid = (TH1D*) GetCoctailMinv(kElId);
+
+ h65cocBgElid->Add(fh_mean_bg_minv[kElId]);
+
+ h65coc1Elid->Add(fh_mean_bg_minv[kElId]);
+ h65coc1Elid->Add(fh_mean_bg_minv[kElId], -1.);
+
+ h65cocBgElid->Rebin(nRebin);
+ h65cocBgElid->Scale(1. / bW);
+ h65coc1Elid->Rebin(nRebin);
+ h65coc1Elid->Scale(1. / bW);
+
+ if (setMinMax) {
+ h65cocBgElid->SetMinimum(yMin);
+ h65cocBgElid->SetMaximum(yMax);
+ }
+
+ h65cocBgElid->GetXaxis()->SetRangeUser(0, 2.);
+ h65cocBgElid->GetYaxis()->SetTitle(yTitle.c_str());
+ h65cocBgElid->SetTitle(0);
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_bgCoc_elid", "minv_CB_6_signal_bgCoc_elid", 800, 800);
+ DrawH1({h65cocBgElid, h65cbAssElid, h65SigCocBgElid, h65coc1Elid}, {"", "", "", ""}, kLinear, kLog, false, 0.8, 0.8,
+ 0.99, 0.99, "hist p");
+
+ TLegend* legend65 = new TLegend(0.55, 0.75, 0.98, 0.9);
+ legend65->SetFillColor(kWhite);
+ legend65->AddEntry(h65cocBgElid, "Cocktail + BG (Elid)");
+ legend65->AddEntry(h65cbAssElid, "B_{c} (assembled, Elid)");
+ legend65->AddEntry(h65SigCocBgElid, "Signal (Coc + BG - B_{c}, Elid)");
+ legend65->AddEntry(h65coc1Elid, "Cocktail (Elid)");
+ legend65->Draw();
+
+ /*fHM[0]->CreateCanvas("minv_CB_6_signal_bgCoc_elid2", "minv_CB_6_signal_bgCoc_elid2", 800, 800);
+ DrawH1({h65SigCocBgElid}, {""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "same pe1");
+
+ TLegend* legend65b = new TLegend(0.55, 0.75, 0.98, 0.9);
+ legend65b->SetFillColor(kWhite);
+ legend65b->AddEntry(h65SigCocBgElid, "Signal (Coc + BG - B_{c}, Elid)");
+ legend65b->Draw();*/
+
+ // Pt cut
+ TH1D* h66SigCocBgPt = (TH1D*) fh_mean_combSignalBCoc_assemb_minv[kPtCut]->Clone();
+ TH1D* h66cbAssPt = (TH1D*) fh_mean_combBg_assemb_minv[kPtCut]->Clone();
+ TH1D* h66cocBgPt = (TH1D*) GetCoctailMinv(kPtCut);
+ h66cocBgPt->Add(fh_mean_bg_minv[kPtCut]);
+
+ h66cocBgPt->Rebin(nRebin);
+ h66cocBgPt->Scale(1. / bW);
+
+ if (setMinMax) {
+ h66cocBgPt->SetMinimum(yMin);
+ h66cocBgPt->SetMaximum(yMax);
+ }
+
+ h66cocBgPt->GetXaxis()->SetRangeUser(0, 2.);
+ h66cocBgPt->GetYaxis()->SetTitle(yTitle.c_str());
+ h66cocBgPt->SetTitle(0);
+
+ fHM[0]->CreateCanvas("minv_CB_6_signal_bgCoc_pt", "minv_CB_6_signal_bgCoc_pt", 800, 800);
+ DrawH1({h66cocBgPt, h66cbAssPt, h66SigCocBgPt}, {"", "", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legend66 = new TLegend(0.55, 0.75, 0.98, 0.9);
+ legend66->SetFillColor(kWhite);
+ legend66->AddEntry(h66cocBgPt, "Cocktail + BG (P_{t})");
+ legend66->AddEntry(h66cbAssPt, "B_{c} (assembled, P_{t})");
+ legend66->AddEntry(h66SigCocBgPt, "Signal (Coc + BG - B_{c}, P_{t})");
+ legend66->Draw();
+
+ /* 7) Error */
+
+ TH1D* h71NpmElid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv_raw[kElId]->Clone();
+ TH1D* h71sigElid = (TH1D*) fh_mean_combSignal_errProp_minv[kElId]->Clone();
+
+ h71NpmElid->GetXaxis()->SetRangeUser(0, 2.);
+ h71NpmElid->GetXaxis()->SetTitle(xTitle.c_str());
+ h71NpmElid->GetYaxis()->SetTitle(yTitle.c_str());
+ h71NpmElid->SetTitle(0);
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double errNpm = h71NpmElid->GetBinError(iBin);
+ h71NpmElid->SetBinContent(iBin, errNpm);
+ }
+
+ fHM[0]->CreateCanvas("minv_CB_7_error", "minv_CB_7_error", 800, 800);
+ DrawH1({h71NpmElid, h71sigElid}, {"", ""}, kLinear, kLog, false, 0.8, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legend71 = new TLegend(0.55, 0.75, 0.98, 0.9);
+ legend71->SetFillColor(kWhite);
+ legend71->AddEntry(h71NpmElid, "error N^{+-}_{same} (Elid)");
+ legend71->AddEntry(h71sigElid, "error Signal = N^{+-}_{same} - B_{c} (Elid)");
+ legend71->Draw();
}
void CbmAnaDielectronTaskDrawAll::DrawSBgVsMinv()
{
TH1D* bg = (TH1D*) fh_mean_bg_minv[kTtCut]->Clone();
+ TH1D* combBg = (TH1D*) fh_mean_combBg_minv[kTtCut]->Clone();
TH1D* coctail = GetCoctailMinv(kTtCut);
+
fh_mean_sbg_vs_minv[kTtCut] =
new TH1D(("fh_sbg_vs_minv_" + CbmLmvmHist::fAnaSteps[kTtCut]).c_str(),
("fh_sbg_vs_minv_" + CbmLmvmHist::fAnaSteps[kTtCut] + ";M_{ee} [GeV/c^{2}];Cocktail/Background").c_str(),
@@ -352,8 +1389,19 @@ void CbmAnaDielectronTaskDrawAll::DrawSBgVsMinv()
fh_mean_sbg_vs_minv[kTtCut]->Scale(1. / 20.);
fh_mean_sbg_vs_minv[kTtCut]->GetXaxis()->SetRangeUser(0, 2.);
+ fh_mean_combSBg_vs_minv[kAcc] =
+ new TH1D(("fh_combSBg_vs_minv_" + CbmLmvmHist::fAnaSteps[kAcc]).c_str(),
+ ("fh_combSBg_vs_minv_" + CbmLmvmHist::fAnaSteps[kAcc] + ";M_{ee} [GeV/c^{2}];Cocktail/comb. BG").c_str(),
+ bg->GetNbinsX(), bg->GetXaxis()->GetXmin(), bg->GetXaxis()->GetXmax());
+ fh_mean_combSBg_vs_minv[kAcc]->Divide(coctail, combBg, 1., 1., "B");
+ fh_mean_combSBg_vs_minv[kAcc]->Rebin(20);
+ fh_mean_combSBg_vs_minv[kAcc]->Scale(1. / 20.);
+ fh_mean_combSBg_vs_minv[kAcc]->GetXaxis()->SetRangeUser(0, 2.);
+
bg = (TH1D*) fh_mean_bg_minv[kPtCut]->Clone();
+ combBg = (TH1D*) fh_mean_combBg_minv[kPtCut]->Clone();
coctail = GetCoctailMinv(kPtCut);
+
fh_mean_sbg_vs_minv[kPtCut] =
new TH1D(("fh_sbg_vs_minv_" + CbmLmvmHist::fAnaSteps[kPtCut]).c_str(),
("fh_sbg_vs_minv_" + CbmLmvmHist::fAnaSteps[kPtCut] + ";M_{ee} [GeV/c^{2}];Cocktail/Background").c_str(),
@@ -363,10 +1411,24 @@ void CbmAnaDielectronTaskDrawAll::DrawSBgVsMinv()
fh_mean_sbg_vs_minv[kPtCut]->Scale(1. / 20.);
fh_mean_sbg_vs_minv[kPtCut]->GetXaxis()->SetRangeUser(0, 2.);
+ /*fh_mean_combSBg_vs_minv[kPtCut] = new TH1D(
+ ("fh_combSBg_vs_minv_" + CbmLmvmHist::fAnaSteps[kPtCut]).c_str(),
+ ("fh_combSBg_vs_minv_" + CbmLmvmHist::fAnaSteps[kPtCut] + ";M_{ee} [GeV/c^{2}];Cocktail/comb. Background").c_str(),
+ bg->GetNbinsX(), bg->GetXaxis()->GetXmin(), bg->GetXaxis()->GetXmax());
+ fh_mean_combSBg_vs_minv[kPtCut]->Divide(coctail, combBg, 1., 1., "B");
+ fh_mean_combSBg_vs_minv[kPtCut]->Rebin(20);
+ fh_mean_combSBg_vs_minv[kPtCut]->Scale(1. / 20.);
+ fh_mean_combSBg_vs_minv[kPtCut]->GetXaxis()->SetRangeUser(0, 2.);*/
+
fHM[0]->CreateCanvas("lmvm_sbg_vs_minv", "lmvm_sbg_vs_minv", 800, 800);
DrawH1({fh_mean_sbg_vs_minv[kTtCut], fh_mean_sbg_vs_minv[kPtCut]}, {"Without Pt cut", "With Pt cut"}, kLinear, kLog,
true, 0.6, 0.85, 0.99, 0.99);
gPad->SetLogy(true);
+
+ /*fHM[0]->CreateCanvas("lmvm_sbgComb_vs_minv", "lmvm_sbgComb_vs_minv", 800, 800);
+ DrawH1({fh_mean_combSBg_vs_minv[kAcc], fh_mean_combSBg_vs_minv[kPtCut]}, {"Without Pt cut", "With Pt cut"}, kLinear,
+ kLog, true, 0.6, 0.85, 0.99, 0.99);*/
+ gPad->SetLogy(true);
}
void CbmAnaDielectronTaskDrawAll::DrawMinvPtAll()
@@ -405,6 +1467,26 @@ void CbmAnaDielectronTaskDrawAll::FillMeanHist()
fh_mean_pi0_minv[step] = (TH1D*) H1(iS, "fh_pi0_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone();
fh_mean_eta_minv_pt[step] = (TH2D*) H2(iS, "fh_eta_minv_pt_" + CbmLmvmHist::fAnaSteps[step])->Clone();
fh_mean_pi0_minv_pt[step] = (TH2D*) H2(iS, "fh_pi0_minv_pt_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsPM_sameEvent_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsPM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsPP_sameEvent_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsPP_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsMM_sameEvent_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsMM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsPM_mixedEvents_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsPM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsPP_mixedEvents_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsPP_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_combPairsMM_mixedEvents_minv[step] =
+ (TH1D*) H1(iS, "fh_combPairsMM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_nof_plutoElectrons[step] =
+ (TH1D*) H1(iS, "fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_nof_plutoPositrons[step] =
+ (TH1D*) H1(iS, "fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_nof_urqmdElectrons[step] =
+ (TH1D*) H1(iS, "fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ fh_mean_nof_urqmdPositrons[step] =
+ (TH1D*) H1(iS, "fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[step])->Clone();
}
else {
fh_mean_bg_minv[step]->Add((TH1D*) H1(iS, "fh_bg_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone());
@@ -412,15 +1494,55 @@ void CbmAnaDielectronTaskDrawAll::FillMeanHist()
fh_mean_pi0_minv[step]->Add((TH1D*) H1(iS, "fh_pi0_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone());
fh_mean_eta_minv_pt[step]->Add((TH2D*) H2(iS, "fh_eta_minv_pt_" + CbmLmvmHist::fAnaSteps[step])->Clone());
fh_mean_pi0_minv_pt[step]->Add((TH2D*) H2(iS, "fh_pi0_minv_pt_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsPM_sameEvent_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsPM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsPP_sameEvent_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsPP_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsMM_sameEvent_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsMM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsPM_mixedEvents_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsPM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsPP_mixedEvents_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsPP_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_combPairsMM_mixedEvents_minv[step]->Add(
+ (TH1D*) H1(iS, "fh_combPairsMM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_nof_plutoElectrons[step]->Add(
+ (TH1D*) H1(iS, "fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_nof_plutoPositrons[step]->Add(
+ (TH1D*) H1(iS, "fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_nof_urqmdElectrons[step]->Add(
+ (TH1D*) H1(iS, "fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[step])->Clone());
+ fh_mean_nof_urqmdPositrons[step]->Add(
+ (TH1D*) H1(iS, "fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[step])->Clone());
}
}
- //TODO: nofSignals without QGP
- fh_mean_bg_minv[step]->Scale(1. / (double) fNofSignals);
- fh_mean_eta_minv[step]->Scale(1. / (double) fNofSignals);
- fh_mean_pi0_minv[step]->Scale(1. / (double) fNofSignals);
- fh_mean_eta_minv_pt[step]->Scale(1. / (double) fNofSignals);
- fh_mean_pi0_minv_pt[step]->Scale(1. / (double) fNofSignals);
- }
+
+ fh_mean_combPairsPM_sameEvent_minv_raw[step] = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[step]->Clone();
+ fh_mean_combPairsPP_sameEvent_minv_raw[step] = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[step]->Clone();
+ fh_mean_combPairsMM_sameEvent_minv_raw[step] = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[step]->Clone();
+ fh_mean_combPairsPM_mixedEvents_minv_raw[step] = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv[step]->Clone();
+ fh_mean_combPairsPP_mixedEvents_minv_raw[step] = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv[step]->Clone();
+ fh_mean_combPairsMM_mixedEvents_minv_raw[step] = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv[step]->Clone();
+
+ Int_t nofEvents = 0;
+ for (int i = 0; i < fNofSignals; i++) {
+ if (!fDrawQgp && i == kQgp) continue;
+ nofEvents += (int) H1(i, "fh_event_number")->GetEntries();
+ cout << "FillMeanHist: nofEvents = " << nofEvents << endl;
+ }
+ cout << "FillMeanHist: final nofEvents = " << nofEvents << endl;
+
+ // Scaling; Comb. BG histograms are scaled in CalcCombBGHistos()
+ fh_mean_bg_minv[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_eta_minv[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_pi0_minv[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_eta_minv_pt[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_pi0_minv_pt[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_nof_plutoElectrons[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_nof_plutoPositrons[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_nof_urqmdElectrons[step]->Scale(1. / (double) (nofEvents));
+ fh_mean_nof_urqmdPositrons[step]->Scale(1. / (double) (nofEvents));
+ } // steps
}
void CbmAnaDielectronTaskDrawAll::SaveHist()
@@ -430,8 +1552,11 @@ void CbmAnaDielectronTaskDrawAll::SaveHist()
TFile* f = TFile::Open(string(fOutputDir + "/draw_all_hist.root").c_str(), "RECREATE");
for (int i = 0; i < CbmLmvmHist::fNofAnaSteps; i++) {
fh_mean_bg_minv[i]->Write();
+ fh_mean_combBg_minv[i]->Write();
+ fh_mean_combBg_GeomMeanSame_minv[i]->Write();
+ fh_mean_combBg_GeomMeanMixed_minv[i]->Write();
+ //fh_mean_combSignalNpm_minv[i]->Write();
fh_mean_eta_minv[i]->Write();
- fh_mean_pi0_minv[i]->Write();
}
fh_mean_sbg_vs_minv[kTtCut]->Write();
fh_mean_sbg_vs_minv[kPtCut]->Write();
@@ -494,6 +1619,293 @@ void CbmAnaDielectronTaskDrawAll::CalcCutEffRange(Double_t minMinv, Double_t max
t->Draw();
}
+void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
+{
+
+ int nofBinsRaw;
+ double bW = 0;
+ Int_t nofEvents = 0;
+ for (int i = 0; i < fNofSignals; i++) {
+ if (!fDrawQgp && i == kQgp) continue;
+ nofEvents += (int) H1(i, "fh_event_number")->GetEntries();
+ cout << "CalcCombBGHistos: nofEvents = " << nofEvents << endl;
+ }
+
+ // first convert MeV into Bin for later normalization
+ int nBinsOrig = 4000; // was set in CbmAnaDielectronTask.cxx
+ int upperLimit = 4000; // in MeV (assumed, that lower limit of histo = 0)
+
+ int rFrom = 300; // lower and upper range to normalize (in Mev/c^2)
+ int rTo = 1000;
+ double binFrom = (double) (1. * nBinsOrig / nRebin) * (1. * rFrom / upperLimit);
+ double binTo = (double) (1. * nBinsOrig / nRebin) * (1. * rTo / upperLimit);
+
+ for (Int_t step = 0; step < CbmLmvmHist::fNofAnaSteps; step++) {
+ fh_mean_combPairsPM_sameEvent_minv[step]->Rebin(nRebin);
+ fh_mean_combPairsPP_sameEvent_minv[step]->Rebin(nRebin);
+ fh_mean_combPairsMM_sameEvent_minv[step]->Rebin(nRebin);
+ fh_mean_combPairsPM_mixedEvents_minv[step]->Rebin(nRebin);
+ fh_mean_combPairsPP_mixedEvents_minv[step]->Rebin(nRebin);
+ fh_mean_combPairsMM_mixedEvents_minv[step]->Rebin(nRebin);
+
+ fh_mean_combPairsPM_sameEvent_minv_raw[step]->Rebin(nRebin);
+ fh_mean_combPairsPP_sameEvent_minv_raw[step]->Rebin(nRebin);
+ fh_mean_combPairsMM_sameEvent_minv_raw[step]->Rebin(nRebin);
+ fh_mean_combPairsPM_mixedEvents_minv_raw[step]->Rebin(nRebin);
+ fh_mean_combPairsPP_mixedEvents_minv_raw[step]->Rebin(nRebin);
+ fh_mean_combPairsMM_mixedEvents_minv_raw[step]->Rebin(nRebin);
+
+ nofBinsRaw = fh_mean_combPairsPM_sameEvent_minv_raw[step]->GetNbinsX();
+
+ // calculate geom. mean of same events
+ TH1D* hBpp = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[step]->Clone();
+ TH1D* hBmm = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[step]->Clone();
+
+ fh_mean_combBg_GeomMeanSame_minv[step] = (TH1D*) hBpp->Clone();
+
+ int nofBins = hBpp->GetNbinsX(); // MIND: 'nofBins' is used throughout this method!
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double m1 = hBpp->GetBinContent(iBin);
+ double m2 = hBmm->GetBinContent(iBin);
+ double content = TMath::Sqrt(m1 * m2);
+ fh_mean_combBg_GeomMeanSame_minv[step]->SetBinContent(iBin, content);
+ }
+
+ // calculate geom. mean of mixed events
+ TH1D* hbpp = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv[step]->Clone();
+ TH1D* hbmm = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv[step]->Clone();
+
+ fh_mean_combBg_GeomMeanMixed_minv[step] = (TH1D*) hbpp->Clone();
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) {
+ double m1 = hbpp->GetBinContent(iBin);
+ double m2 = hbmm->GetBinContent(iBin);
+ double content = 0;
+ if (m1 == 0 && m2 != 0) content = m2;
+ else if (m1 != 0 && m2 == 0)
+ content = m1;
+ else
+ content = TMath::Sqrt(m1 * m2);
+ fh_mean_combBg_GeomMeanMixed_minv[step]->SetBinContent(iBin, content);
+ }
+
+ // normalization factor for same/mixed geom. mean
+ double intSameGM = fh_mean_combBg_GeomMeanSame_minv[step]->Integral(binFrom, binTo);
+ double intMixedGM = fh_mean_combBg_GeomMeanMixed_minv[step]->Integral(binFrom, binTo);
+ double normGM = (double) intSameGM / intMixedGM;
+ cout << "step " << step << ": normGM = " << normGM << endl;
+
+ // calculate k factor
+ TH1D* k = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv[step]->Clone();
+ k->Divide(fh_mean_combBg_GeomMeanMixed_minv[step]);
+ k->Scale(1. / 2);
+ fh_mean_combBg_k_minv[step] = (TH1D*) k->Clone();
+
+ // calculate combinatorial BG
+ TH1D* Bc = (TH1D*) k->Clone();
+ Bc->Multiply(fh_mean_combBg_GeomMeanSame_minv[step]);
+ Bc->Scale(2.);
+ fh_mean_combBg_minv[step] = (TH1D*) Bc->Clone();
+ fh_mean_combBg_raw_minv[step] = (TH1D*) Bc->Clone();
+ fh_mean_combBg_raw_minv[step]->Scale(fNofSignals * nofEvents);
+
+ // calculate assembled combinatorial BG from same (<= 0.3 GeV) and mixed (> 0.3 GeV) events data
+ fh_mean_combBg_assemb_minv[step] = (TH1D*) fh_mean_combBg_minv[step]->Clone();
+
+ for (int iBin = binFrom + 1; iBin <= nofBins; iBin++) { // from > 300 MeV on normalized data from mixed event
+ double m1 = fh_mean_combBg_k_minv[step]->GetBinContent(iBin);
+ double m2 = fh_mean_combBg_GeomMeanMixed_minv[step]->GetBinContent(iBin);
+ double content = 2 * m1 * m2 * normGM;
+ fh_mean_combBg_assemb_minv[step]->SetBinContent(iBin, content);
+ }
+
+ // calculate comb. signal
+ // from 'N+-same'
+ fh_mean_combSignalNpm_minv[step] = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[step]->Clone();
+ fh_mean_combSignalNpm_minv[step]->Add(fh_mean_combBg_minv[step], -1.);
+
+ TH1D* cock;
+
+ // from 'Cocktail + BG'
+ if (step == kMc) cock = (TH1D*) GetCoctailMinv(kMc);
+ else if (step == kAcc)
+ cock = (TH1D*) GetCoctailMinv(kAcc);
+ else if (step == kReco)
+ cock = (TH1D*) GetCoctailMinv(kReco);
+ else if (step == kChi2Prim)
+ cock = (TH1D*) GetCoctailMinv(kChi2Prim);
+ else if (step == kElId)
+ cock = (TH1D*) GetCoctailMinv(kElId);
+ else if (step == kGammaCut)
+ cock = (TH1D*) GetCoctailMinv(kGammaCut);
+ else if (step == kMvd1Cut)
+ cock = (TH1D*) GetCoctailMinv(kMvd1Cut);
+ else if (step == kMvd2Cut)
+ cock = (TH1D*) GetCoctailMinv(kMvd2Cut);
+ else if (step == kStCut)
+ cock = (TH1D*) GetCoctailMinv(kStCut);
+ else if (step == kRtCut)
+ cock = (TH1D*) GetCoctailMinv(kRtCut);
+ else if (step == kTtCut)
+ cock = (TH1D*) GetCoctailMinv(kTtCut);
+ else if (step == kPtCut)
+ cock = (TH1D*) GetCoctailMinv(kPtCut);
+
+ cock->Rebin(nRebin);
+ fh_mean_combSignalNpm_minv[step] = (TH1D*) fh_mean_bg_minv[step]->Clone();
+ fh_mean_combSignalNpm_minv[step]->Add(cock);
+ fh_mean_combSignalNpm_minv[step]->Add(fh_mean_combBg_minv[step], -1.);
+
+ // calculate assembled comb. signal from same (<= 0.3 GeV) and mixed (> 0.3 GeV) events data
+ // from 'N+-same'
+ fh_mean_combSignalNpm_assemb_minv[step] = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[step]->Clone();
+ fh_mean_combSignalNpm_assemb_minv[step]->Add(fh_mean_combBg_assemb_minv[step], -1.);
+ fh_mean_combSignalNpm_assemb_minv[step]->SetTitle("Signal (N^{+-}_{same} - B_{C})");
+
+ // from 'Cocktail + BG'
+ fh_mean_combSignalBCoc_assemb_minv[step] = (TH1D*) fh_mean_bg_minv[step]->Clone();
+ fh_mean_combSignalBCoc_assemb_minv[step]->Add(cock);
+ TH1D* cbAss = (TH1D*) fh_mean_combBg_assemb_minv[step]->Clone();
+ cbAss->Scale(1. / nofEvents);
+ fh_mean_combSignalBCoc_assemb_minv[step]->Add(cbAss, -1.);
+
+ // error calculation
+ TH1D* err = (TH1D*) fh_mean_combPairsPM_sameEvent_minv_raw[step]->Clone();
+ err->Reset("ice");
+ fh_mean_combBg_errProp_minv[step] = (TH1D*) err->Clone();
+ fh_mean_combSignal_errProp_minv[step] = (TH1D*) err->Clone();
+
+ for (int iBin = 1; iBin <= nofBinsRaw; iBin++) { // calculate error propagation
+ double Npm = fh_mean_combPairsPM_sameEvent_minv_raw[step]->GetBinContent(iBin);
+ double Bpp = fh_mean_combPairsPP_sameEvent_minv_raw[step]->GetBinContent(iBin);
+ double Bmm = fh_mean_combPairsMM_sameEvent_minv_raw[step]->GetBinContent(iBin);
+ double bpm = fh_mean_combPairsPM_mixedEvents_minv_raw[step]->GetBinContent(iBin);
+ double bpp = fh_mean_combPairsPP_mixedEvents_minv_raw[step]->GetBinContent(iBin);
+ double bmm = fh_mean_combPairsMM_mixedEvents_minv_raw[step]->GetBinContent(iBin);
+
+ // calculation of error propagation
+ double DNpm = TMath::Sqrt(Npm); // Δ<B+->
+ double DBpp = TMath::Sqrt(Bpp); // Δ<B++>
+ double DBmm = TMath::Sqrt(Bmm); // Δ<B-->
+ double Dbpm = TMath::Sqrt(bpm); // Δ<b+->
+ double Dbpp = TMath::Sqrt(bpp); // Δ<b++>
+ double Dbmm = TMath::Sqrt(bmm); // Δ<b-->
+
+ double dNpm = (double) 1.; // derivatives of B_c and signal resp. to single contributions
+ double dBpp = (double) (1. / 2) * bpm * Bmm * (1. / TMath::Sqrt(Bpp * Bmm * bpp * bmm));
+ double dBmm = (double) (1. / 2) * bpm * Bpp * (1. / TMath::Sqrt(Bpp * Bmm * bpp * bmm));
+ double dbpm = (double) TMath::Sqrt((Bpp * Bmm) / (bpp * bmm));
+ double dbpp = (double) ((-1.) / 2) * bpm * bmm * TMath::Sqrt(Bpp * Bmm)
+ / (TMath::Sqrt(bpp * bmm) * TMath::Sqrt(bpp * bmm) * TMath::Sqrt(bpp * bmm));
+ double dbmm = (double) ((-1.) / 2) * bpm * bpp * TMath::Sqrt(Bpp * Bmm)
+ / (TMath::Sqrt(bpp * bmm) * TMath::Sqrt(bpp * bmm) * TMath::Sqrt(bpp * bmm));
+ double dbpm2 = (double) normGM; // for > 300 MeV
+
+ double fNpm = (DNpm * dNpm) * (DNpm * dNpm);
+ double fBpp = (DBpp * dBpp) * (DBpp * dBpp); // single contribution factors to error propagation
+ double fBmm = (DBmm * dBmm) * (DBmm * dBmm);
+ double fbpm = (Dbpm * dbpm) * (Dbpm * dbpm);
+ double fbpp = (Dbpp * dbpp) * (Dbpp * dbpp); // for > 300 MeV
+ double fbmm = (Dbmm * dbmm) * (Dbmm * dbmm);
+ double fbpm2 = (Dbpm * dbpm2) * (Dbpm * dbpm2); // for > 300 MeV
+
+ double errorBc = TMath::Sqrt(fBpp + fBmm + fbpm + fbpp + fbmm); // final error propagation value
+ double errorBc2 = TMath::Sqrt(fbpm2);
+
+ double errorSig = TMath::Sqrt(fNpm + fBpp + fBmm + fbpm + fbpp + fbmm);
+ double errorSig2 = TMath::Sqrt(fNpm + fbpm2);
+
+ // fill error histograms (histos with error information only)
+ if (iBin <= binFrom) fh_mean_combBg_errProp_minv[step]->SetBinContent(iBin, errorBc);
+ if (iBin > binFrom) fh_mean_combBg_errProp_minv[step]->SetBinContent(iBin, errorBc2);
+ if (iBin <= binFrom) fh_mean_combSignal_errProp_minv[step]->SetBinContent(iBin, errorSig);
+ if (iBin > binFrom) fh_mean_combSignal_errProp_minv[step]->SetBinContent(iBin, errorSig2);
+
+ // set error value in CB histograms
+ fh_mean_combBg_minv[step]->SetBinError(iBin, errorBc);
+ if (iBin <= binFrom) fh_mean_combBg_assemb_minv[step]->SetBinError(iBin, errorBc);
+ if (iBin > binFrom) fh_mean_combBg_assemb_minv[step]->SetBinError(iBin, errorBc2);
+ if (iBin <= binFrom) fh_mean_combSignalNpm_assemb_minv[step]->SetBinError(iBin, errorSig);
+ if (iBin > binFrom) fh_mean_combSignalNpm_assemb_minv[step]->SetBinError(iBin, errorSig2);
+
+ if (iBin % (500 / nRebin) == 0 && iBin <= (2000 / nRebin)) {
+ cout << "step = " << step << endl;
+ cout << "Npm = " << Npm << endl;
+ cout << "bpm = " << bpm << endl;
+ cout << "bp = " << bpp << endl;
+ cout << "bm = " << bmm << endl;
+ cout << "Bp = " << Bpp << endl;
+ cout << "Bm = " << Bmm << endl;
+ //cout << "k = " << k << endl;
+ cout << "Dbpm = " << Dbpm << endl;
+ cout << "Dbp = " << Dbpp << endl;
+ cout << "Dbm = " << Dbmm << endl;
+ cout << "DBp = " << DBpp << endl;
+ cout << "DBm = " << DBmm << endl;
+ cout << "dbpm = " << dbpm << endl;
+ cout << "dbp = " << dbpp << endl;
+ cout << "dbm = " << dbmm << endl;
+ cout << "dBp = " << dBpp << endl;
+ cout << "dBm = " << dBmm << endl;
+ cout << "fbpm = " << fbpm << endl;
+ cout << "fbp = " << fbpp << endl;
+ cout << "fbm = " << fbmm << endl;
+ cout << "fBp = " << fBpp << endl;
+ cout << "fBm = " << fBmm << endl;
+ cout << "errPropBc = " << errorBc << endl;
+ cout << "errPropBc2 = " << errorBc2 << endl;
+ cout << "errorSig = " << errorSig << endl;
+ cout << "errorSig2 = " << errorSig2 << endl;
+ cout << "Bin error CB = " << fh_mean_combBg_errProp_minv[step]->GetBinContent(iBin) << endl;
+ cout << "Bin error Signal = " << fh_mean_combSignal_errProp_minv[step]->GetBinContent(iBin) << endl;
+ cout << "Bin error Npm = " << fh_mean_combPairsPM_sameEvent_minv_raw[step]->GetBinError(iBin) << endl;
+ }
+ } // error propagation
+
+ // scale histograms
+ bW = fh_mean_combPairsPM_sameEvent_minv[step]->GetBinWidth(1);
+ fh_mean_combPairsPM_sameEvent_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combPairsPP_sameEvent_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combPairsMM_sameEvent_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combPairsPM_mixedEvents_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combPairsPP_mixedEvents_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combPairsMM_mixedEvents_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combBg_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combBg_assemb_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combBg_GeomMeanSame_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combBg_GeomMeanMixed_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combSignalNpm_minv[step]->Scale(1. / (nofEvents * bW));
+ fh_mean_combSignalNpm_assemb_minv[step]->Scale(1. / (nofEvents * bW));
+ //fh_mean_combSignalBCoc_minv[step]->Scale(1. / (nofEvents * bW)); // is this histo needed?
+ fh_mean_combSignalBCoc_assemb_minv[step]->Scale(1. / (bW)); // had been normalized earlier due to Cocktail
+ } // steps
+}
+
+/*void CbmAnaDielectronTaskDrawAll::CompareSTSversions()
+{
+ TFile* _Input_file_v16g = TFile::Open("/lustre/cbm/users/criesen/cbm/data/lmvm/inmed/analysis.all.root");
+ TFile* _Input_file_v19a = TFile::Open("/lustre/cbm/users/criesen/cbm/data/lmvm_sts-v19a_5M/inmed/analysis.all.root");
+
+ TH1D* hv16g = (TH1D*) _Input_file_v16g->Get("fh_signal_minv_mc");
+ TH1D* hv19a = (TH1D*) _Input_file_v19a->Get("fh_signal_minv_mc");
+
+ hv16g->Divide(hv19a);
+
+ hv16g->Rebin(nRebin);
+ hv16g->Scale(1. / nRebin);
+
+ hv16g->GetXaxis()->SetRangeUser(0, 2.);
+ hv16g->GetXaxis()->SetTitle("M_{ee} [GeV/c^{2}]");
+
+ hv16g->SetTitle(0);
+
+ DrawH1({hv16g}, {""}, kLinear, kLog, false, 0.8, 0.8, 0.95, 0.95, "HIST L");
+ hv16g->GetYaxis()->SetTitle("ratio v16g/v19a");
+ hv16g->GetYaxis()->SetLabelSize(0.05);
+
+ hv16g->Draw();
+}*/
TH1D* CbmAnaDielectronTaskDrawAll::SBgRange(Double_t min, Double_t max)
{
@@ -578,8 +1990,10 @@ void CbmAnaDielectronTaskDrawAll::DrawSBgSignals()
if (step < kElId) continue;
if (!fUseMvd && (step == kMvd1Cut || step == kMvd2Cut)) { continue; }
- TH1D* s = (TH1D*) H1(iF, "fh_signal_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone();
- TH1D* bg = (TH1D*) fh_mean_bg_minv[step]->Clone();
+ TH1D* s = (TH1D*) H1(iF, "fh_signal_minv_" + CbmLmvmHist::fAnaSteps[step])->Clone();
+ TH1D* bg = (TH1D*) fh_mean_bg_minv[step]->Clone();
+ Int_t nofEvents = (int) H1(iF, "fh_event_number")->GetEntries();
+ s->Scale(1. / nofEvents);
cFit->cd();
if (iF == kPhi) {
if (s->GetEntries() > 0) s->Fit("gaus", "Q", "", 0.95, 1.05);
@@ -637,7 +2051,7 @@ void CbmAnaDielectronTaskDrawAll::DrawSBgSignals()
void CbmAnaDielectronTaskDrawAll::SaveCanvasToImage()
{
cout << "Images output dir:" << fOutputDir << endl;
- fHM[0]->SaveCanvasToImage(fOutputDir, "eps;png");
+ fHM[0]->SaveCanvasToImage(fOutputDir, "png"); // fHM[0]->SaveCanvasToImage(fOutputDir, "eps;png");
}
ClassImp(CbmAnaDielectronTaskDrawAll);
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.h b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.h
old mode 100644
new mode 100755
index de0f9c7e928db814ff8a7cd3b7225d898e3e4d81..6476b6b51d05419a14a108fc999d3ee8101222fc
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.h
@@ -81,6 +81,8 @@ public:
private:
static const int fNofSignals = 5;
+ int nRebin;
+
Bool_t fUseMvd; // do you want to draw histograms related to the MVD detector?
Bool_t fDrawQgp; // do you wan to draq QGP signal
@@ -95,6 +97,41 @@ private:
std::vector<TH2D*> fh_mean_pi0_minv_pt;
std::vector<TH1D*> fh_mean_sbg_vs_minv; //Coctail/BG vs. invariant mass for different analysis steps
+ // Combinatorial histograms
+ std::vector<TH1D*> fh_mean_combPairsPM_sameEvent_minv;
+ std::vector<TH1D*> fh_mean_combPairsPP_sameEvent_minv;
+ std::vector<TH1D*> fh_mean_combPairsMM_sameEvent_minv;
+ std::vector<TH1D*> fh_mean_combPairsPM_mixedEvents_minv;
+ std::vector<TH1D*> fh_mean_combPairsPP_mixedEvents_minv;
+ std::vector<TH1D*> fh_mean_combPairsMM_mixedEvents_minv;
+ std::vector<TH1D*> fh_mean_combPairsPM_sameEvent_minv_raw; // won't be scaled with binWidth
+ std::vector<TH1D*> fh_mean_combPairsPP_sameEvent_minv_raw;
+ std::vector<TH1D*> fh_mean_combPairsMM_sameEvent_minv_raw;
+ std::vector<TH1D*> fh_mean_combPairsPM_mixedEvents_minv_raw;
+ std::vector<TH1D*> fh_mean_combPairsPP_mixedEvents_minv_raw;
+ std::vector<TH1D*> fh_mean_combPairsMM_mixedEvents_minv_raw;
+ std::vector<TH1D*> fh_mean_combBg_errProp_minv; // to observe error propagation
+ std::vector<TH1D*> fh_mean_combSignal_errProp_minv;
+ std::vector<TH1D*> fh_mean_combBg_GeomMeanSame_minv; // geom. mean of comb. BG ( := SQRT[(B++) * (B--)] )
+ std::vector<TH1D*> fh_mean_combBg_GeomMeanMixed_minv; // geom. mean of comb. BG ( := SQRT[(b++) * (b--)] )
+ std::vector<TH1D*> fh_mean_combBg_k_minv; // k = (b+-) / ( 2 * Sqrt[(b++) * (b--)] )
+ std::vector<TH1D*> fh_mean_combBg_minv; // combinatorial BG ( := B = 2 * geomMean * k )
+ std::vector<TH1D*> fh_mean_combBg_raw_minv; // won't be scaled; for error propagation
+ std::vector<TH1D*>
+ fh_mean_combBg_assemb_minv; // as previous, but with assembled same (0 - 0.3 GeV) and mixed (0.3 - 2 GeV) data
+ std::vector<TH1D*> fh_mean_combSignalNpm_minv; // combinatorial signal ( := cSig = (N+-) - B )
+ std::vector<TH1D*>
+ fh_mean_combSignalNpm_assemb_minv; // as previous, but with assembled same (0 - 0.3 GeV) and mixed (0.3 - 2 GeV) data
+ std::vector<TH1D*> fh_mean_combSignalBCoc_minv; // combinatorial signal ( := cSig = (Coc + BG) - B )
+ std::vector<TH1D*>
+ fh_mean_combSignalBCoc_assemb_minv; // as previous, but with assembled same (0 - 0.3 GeV) and mixed (0.3 - 2 GeV) data
+ std::vector<TH1D*> fh_mean_combSBg_vs_minv; // cocktail/combBG
+
+ // Number of charged particles vs. momentum
+ std::vector<TH1D*> fh_mean_nof_plutoElectrons;
+ std::vector<TH1D*> fh_mean_nof_plutoPositrons;
+ std::vector<TH1D*> fh_mean_nof_urqmdElectrons;
+ std::vector<TH1D*> fh_mean_nof_urqmdPositrons;
// index: AnalysisSteps
std::vector<TH1D*> fh_sum_s_minv; // sum of all signals
@@ -126,6 +163,12 @@ private:
*/
void DrawMinv(CbmLmvmAnalysisSteps step);
+ /**
+ * \brief Draw invariant mass spectra for all signal types for specified analysis step with BG reduced by combinatorial BG.
+ * \param[in] step Analysis step.
+ */
+ void DrawMinvCombSignalAndBg();
+
/**
* \brief Draw invariant mass vs Pt histograms.
*/
@@ -159,6 +202,16 @@ private:
*/
void CalcCutEffRange(Double_t minMinv, Double_t maxMinv);
+ /**
+ * \brief Calculate combinatorial BG contribution.
+ */
+ void CalcCombBGHistos();
+
+ /**
+ * \brief To compare outputs from simulations with different STS versions
+ */
+ void CompareSTSversions();
+
/**
* \brief Create S/BG vs cuts for specified invariant mass range.
* \param[in] min Minimum invariant mass.
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCandidate.h b/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCandidate.h
index 60e9ff1399b0b856c63486c864be998ed8e8cd2b..144725dc3f1b4c10c5c3ccc9a5ac247ef8495677 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCandidate.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCandidate.h
@@ -24,6 +24,7 @@ public:
, fChi2Prim(0.)
, fChi2sts(0.)
, fMcMotherId(-1)
+ , fEventNumber(0.)
, fStsMcTrackId(-1)
, fRichMcTrackId(-1)
, fTrdMcTrackId(-1)
@@ -75,6 +76,7 @@ public:
Double_t fChi2sts;
Int_t fMcMotherId;
+ Int_t fEventNumber;
Int_t fStsMcTrackId;
Int_t fRichMcTrackId;
Int_t fTrdMcTrackId;
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCuts.h b/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCuts.h
index 3b7b28e9e1c090bb1eb0422f0866f1412ea39bd9..9b7497116cc2ad9025b6c468a4a88f1ab9a9c2b1 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCuts.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmLmvmCuts.h
@@ -57,8 +57,8 @@ public:
fStCutPP = 1.;
//fTtCutAngle = 0.75;
//fTtCutPP = 4.0;
- fTtCutAngle = 1.5;
- fTtCutPP = 1.7;
+ fTtCutAngle = 2.0;
+ fTtCutPP = 3.0;
//fRtCutAngle = 1.0;
//fRtCutPP = 2.5;
fRtCutAngle = 1.2;
diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py
old mode 100755
new mode 100644
index 1630a920a3598e5367c5f8b16e02b26e738d160a..f76179266b26a1f0b9f15b45778e351459acecfd
--- a/macro/analysis/dielectron/batch_job.py
+++ b/macro/analysis/dielectron/batch_job.py
@@ -1,76 +1,38 @@
-#!/usr/bin/env python3
+#!/ usr / bin / env python3
-import os # operating system interfaces
-import sys # provides access to some variables used or maintained by the interpreter and to functions that interact strongly with the interpreter
-import shutil # offers a number of high-level operations on files and collections of files
+import os #operating system interfaces import sys #provides access to some variables used or maintained by the interpreter and to functions that interact strongly with the interpreter import shutil #offers a number of high - level operations on files and collections of files
-def main():
- dataDir = sys.argv[1]
- geoSetup = sys.argv[2]
- plutoParticle = sys.argv[3]
- cbmrootConfigPath = "/lustre/nyx/cbm/users/criesen/build/config.sh"
- macroDir = "/lustre/nyx/cbm/users/criesen/cbmroot/macro/analysis/dielectron/"
- nofEvents = 1000
-
- taskId = os.environ.get('SLURM_ARRAY_TASK_ID')
+ def main() : dataDir = sys.argv[1] geoSetup = sys.argv[2] plutoParticle = sys.argv[3] cbmrootConfigPath = "/lustre/nyx/cbm/users/criesen/build/config.sh" macroDir = "/lustre/nyx/cbm/users/criesen/cbmroot/macro/analysis/dielectron/" nofEvents = 1000
- jobId = os.environ.get('SLURM_ARRAY_JOB_ID')
- colEnergy = "8gev"
- colSystem = "auau"
- print("dataDir:" + dataDir)
+ taskId = os.environ.get('SLURM_ARRAY_TASK_ID')
- os.system(("source {}").format(cbmrootConfigPath))
+ jobId = os.environ.get('SLURM_ARRAY_JOB_ID') colEnergy = "8gev" colSystem = "auau" print("dataDir:" + dataDir)
- workDir = dataDir + "/workdir/" + jobId + "_" + taskId + "/"
- if os.path.exists(workDir):
- shutil.rmtree(workDir)
- os.makedirs(workDir)
- os.chdir(workDir)
+ os.system(("source {}").format(cbmrootConfigPath))
- plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId)
- #plutoFile = getPlutoPath("auau", "8gev", plutoParticle, taskId)
-
- urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/"+colEnergy+"/centr/urqmd.auau.8gev.centr." + str(taskId).zfill(5) + ".root"
- #urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/8gev/centr/urqmd.auau.8gev.centr.00001.root"
-
- mcFile = dataDir + "/mc." + taskId + ".root"
- parFile = dataDir + "/param." + taskId + ".root"
- digiFile = dataDir + "/digi." + taskId + ".root"
- recoFile = dataDir + "/reco." + taskId + ".root"
- litQaFile = dataDir + "/litqa." + taskId + ".root"
- analysisFile = dataDir + "/analysis." + taskId + ".root"
- geoSimFile = dataDir + "/geosim." + taskId + ".root"
- os.system(('root -l -b -q {}/run_sim.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
- macroDir, urqmdFile, plutoFile, mcFile, parFile, geoSimFile, geoSetup, nofEvents))
- os.system(('root -l -b -q {}/run_digi.C\(\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
- macroDir, mcFile, parFile, digiFile, nofEvents))
- os.system(('root -l -b -q {}/run_reco.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
- macroDir, mcFile, parFile, digiFile, recoFile, geoSetup, nofEvents))
- os.system(('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
- macroDir, mcFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents))
- os.system(('root -l -b -q {}/run_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
- macroDir, mcFile, parFile, digiFile, recoFile, analysisFile, plutoParticle, colSystem, colEnergy, geoSetup, nofEvents))
+ workDir = dataDir + "/workdir/" + jobId + "_" + taskId + "/" if os.path.exists(workDir) : shutil.rmtree(workDir) os.makedirs(workDir) os.chdir(workDir)
+#plutoFile = ""
+ plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId) #this one was activated before
+#plutoFile = getPlutoPath("auau", "8gev", plutoParticle, taskId)
-def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId):
- if plutoParticle == "rho0":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/"+colEnergy+"/rho0/epem/pluto.auau."+colEnergy+".rho0.epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "omegaepem":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/"+colEnergy+"/omega/epem/pluto.auau."+colEnergy+".omega.epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "omegadalitz":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/"+colEnergy+"/omega/pi0epem/pluto.auau."+colEnergy+".omega.pi0epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "phi":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/"+colEnergy+"/phi/epem/pluto.auau."+colEnergy+".phi.epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "pi0":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/"+colEnergy+"/pi0/gepem/pluto.auau."+colEnergy+".pi0.gepem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "inmed":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/"+colEnergy+"/rapp.inmed/epem/pluto.auau."+colEnergy+".rapp.inmed.epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "qgp":
- return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/"+colEnergy+"/rapp.qgp/epem/pluto.auau."+colEnergy+".rapp.qgp.epem." + str(taskId).zfill(4) + ".root"
- elif plutoParticle == "urqmd":
- return ""
+ urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/" + colEnergy + "/centr/urqmd.auau.8gev.centr." + str(taskId).zfill(5) + ".root"
+#urqmdFile = "/lustre/nyx/cbm/prod/gen/urqmd/auau/8gev/centr/urqmd.auau.8gev.centr.00001.root"
+ mcFile = dataDir + "/mc." + taskId + ".root" parFile = dataDir + "/param." + taskId + ".root" digiFile = dataDir + "/digi." + taskId + ".root" recoFile = dataDir + "/reco." + taskId + ".root" litQaFile = dataDir + "/litqa." + taskId + ".root" analysisFile = dataDir + "/analysis." + taskId + ".root" geoSimFile = dataDir + "/geosim." + taskId + ".root"
+#os.system(('root -l -b -q {}/run_sim.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
+#macroDir, urqmdFile, plutoFile, mcFile, parFile, geoSimFile, geoSetup, nofEvents))
+#os.system(('root -l -b -q {}/run_digi.C\(\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
+#macroDir, mcFile, parFile, digiFile, nofEvents))
+#os.system(('root -l -b -q {}/run_reco.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
+#macroDir, mcFile, parFile, digiFile, recoFile, geoSetup, nofEvents))
+#os.system(('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(
+#macroDir, mcFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents))
+ os.system(('root -l -b -q {}/run_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, analysisFile, plutoParticle, colSystem, colEnergy, geoSetup, nofEvents))
-if __name__ == '__main__':
- main()
+
+ def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId) : if plutoParticle == "rho0" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/rho0/epem/pluto.auau." + colEnergy + ".rho0.epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "omegaepem" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/omega/epem/pluto.auau." + colEnergy + ".omega.epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "omegadalitz" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/omega/pi0epem/pluto.auau." + colEnergy + ".omega.pi0epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "phi" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/phi/epem/pluto.auau." + colEnergy + ".phi.epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "pi0" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/pi0/gepem/pluto.auau." + colEnergy + ".pi0.gepem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "inmed" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/" + colEnergy + "/rapp.inmed/epem/pluto.auau." + colEnergy + ".rapp.inmed.epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "qgp" : return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/" + colEnergy + "/rapp.qgp/epem/pluto.auau." + colEnergy + ".rapp.qgp.epem." + str(taskId).zfill(4) + ".root" elif plutoParticle == "urqmd" : return ""
+
+
+ if __name__ == '__main__' : main()
diff --git a/macro/analysis/dielectron/batch_job_common.py b/macro/analysis/dielectron/batch_job_common.py
new file mode 100644
index 0000000000000000000000000000000000000000..2e49d9406f9dd0a57235de2f75e501cc5958a708
--- /dev/null
+++ b/macro/analysis/dielectron/batch_job_common.py
@@ -0,0 +1,58 @@
+#!/ usr / bin / env python3
+
+import os #operating system interfaces import sys #provides access to some variables used or maintained by the interpreter and to functions that interact strongly with the interpreter import shutil #offers a number of high - level operations on files and collections of files
+
+
+ def main() : recoDataDirIn = sys.argv[1] dataDirOut = sys.argv[2] geoSetup = sys.argv[3] plutoParticle = sys.argv[4]
+
+ cbmrootConfigPath = "/lustre/nyx/cbm/users/criesen/build/config.sh" macroDir = "/lustre/nyx/cbm/users/criesen/cbmroot/macro/analysis/dielectron/" nofEvents = 100
+
+ taskId = os.environ.get('SLURM_ARRAY_TASK_ID')
+
+ jobId = os.environ.get('SLURM_ARRAY_JOB_ID') colEnergy = "12gev" colEnergyFile = "12agev" colSystem = "auau" print("recoDataDirIn:" + recoDataDirIn) print("dataDirOut:" + dataDirOut)
+
+ os.system(("source {}").format(cbmrootConfigPath))
+
+ workDirOut = dataDirOut + "/workdir/" + jobId + "_" + taskId + "/" if os.path.exists(workDirOut) : shutil.rmtree(workDirOut) os.makedirs(workDirOut) os.chdir(workDirOut)
+
+#offset describes the name of the folder with lowest number - 1 in dataIn directory
+ if plutoParticle == "inmed_had_epem" : offset = 5000 offsetId = offset + int(taskId) if plutoParticle == "w" : offset = 0 offsetId = offset + int(taskId) if plutoParticle == "wdalitz" : offset = 1000 offsetId = offset + int(taskId) if plutoParticle == "phi" : offset = 3000 offsetId = offset + int(taskId) if plutoParticle == "qgp_epem" : offset = 7500 offsetId = offset + int(taskId)
+
+#unlike others, paths of w and wdalitz input do not contain '/TGeant4/'
+ if plutoParticle == "w" or plutoParticle == "wdalitz" : mcFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".event.tra.root" parFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".event.par.root" digiFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".event.raw.root" recoFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".event.rec.root" else : mcFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".event.tra.root" parFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".event.par.root" digiFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".event.raw.root" recoFile = recoDataDirIn + colSystem + "/" + colEnergyFile + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".event.rec.root"
+
+ litQaFile = dataDirOut + "/litqa." + taskId + ".root" analysisFile = dataDirOut + "/analysis." + taskId + ".root"
+#geoSimFile = dataDir + "/geosim." + taskId + ".root"##needed ? ##
+
+#it is assumed to take mc, digi, par and reco files from common folder, so only litqa and ana have to be run
+ os.system(('root -l -b -q {}/run_litqa.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, litQaFile, geoSetup, nofEvents)) os.system(('root -l -b -q {}/run_analysis.C\(\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",\\"{}\\",{}\)').format(macroDir, mcFile, parFile, digiFile, recoFile, analysisFile, plutoParticle, colSystem, colEnergy, geoSetup, nofEvents))
+
+
+#def getPlutoPath(colSystem, colEnergy, plutoParticle, taskId):
+#if plutoParticle == "rho0":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/rho0/epem/pluto.auau." + colEnergy + ".rho0.epem." \
+ + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "omegaepem":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/omega/epem/pluto.auau." + colEnergy \
+ + ".omega.epem." + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "omegadalitz":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/omega/pi0epem/pluto.auau." + colEnergy \
+ + ".omega.pi0epem." + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "phi":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/phi/epem/pluto.auau." + colEnergy + ".phi.epem." \
+ + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "pi0":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktA/" + colEnergy + "/pi0/gepem/pluto.auau." + colEnergy + ".pi0.gepem." \
+ + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "inmed":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/" + colEnergy + "/rapp.inmed/epem/pluto.auau." + colEnergy \
+ + ".rapp.inmed.epem." + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "qgp":
+#return "/lustre/nyx/cbm/prod/gen/pluto/auau/cktRapp/" + colEnergy + "/rapp.qgp/epem/pluto.auau." + colEnergy \
+ + ".rapp.qgp.epem." + str(taskId).zfill(4) + ".root"
+#elif plutoParticle == "urqmd":
+#return ""
+
+
+if __name__ == '__main__':
+ main()
diff --git a/macro/analysis/dielectron/batch_send.py b/macro/analysis/dielectron/batch_send.py
old mode 100755
new mode 100644
index d7686adf0809692244555f79541c315be8c793a8..60223be59bae3a4eb99a198002522aefd22754c6
--- a/macro/analysis/dielectron/batch_send.py
+++ b/macro/analysis/dielectron/batch_send.py
@@ -1,54 +1,29 @@
-#!/usr/bin/env python3
-
-import os
-import shutil
-
-def main():
- nofJobs = 10
- timeLimit="08:00:00"
- geoSetup="sis100_electron"
- plutoParticles = ["inmed", "omegadalitz", "phi", "omegaepem", "qgp"];
-
- #ADDED
- mirrorRotation=15
- dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/"
- jobName = "RICH" + str(mirrorRotation)
-
- #All data in data dir will be removed
- removeData = False
- #jobName = "LMVM"
- if removeData:
- print("All data in dataDir will be removed. Dir:" + dataDir)
- print("Removing...")
- if os.path.exists(dataDir):
- shutil.rmtree(dataDir)
- os.makedirs(dataDir,exist_ok=True)
-
- for plutoParticle in plutoParticles:
- dataDirPluto = dataDir + plutoParticle
- logFile = dataDirPluto + "/log/log_slurm-%A_%a.out"
- errorFile = dataDirPluto + "/error/error_slurm-%A_%a.out"
- workDir = dataDirPluto + "/workdir/"
- makeLogErrorDirs(logFile, errorFile, workDir)
-
- #-p debug
- commandStr=('sbatch --job-name={} --time={} --output={} --error={} --array=1-{} batch_job.py {} {} {}').format(jobName, timeLimit, logFile, errorFile, nofJobs, dataDirPluto, geoSetup, plutoParticle)
-
- os.system(commandStr)
-
-
-def makeLogErrorDirs(logFile, errorFile, workDir):
- if os.path.exists(os.path.dirname(logFile)):
- shutil.rmtree(os.path.dirname(logFile))
- os.makedirs(os.path.dirname(logFile),exist_ok=True)
-
- if os.path.exists(os.path.dirname(errorFile)):
- shutil.rmtree(os.path.dirname(errorFile))
- os.makedirs(os.path.dirname(errorFile),exist_ok=True)
-
- if os.path.exists(os.path.dirname(workDir)):
- shutil.rmtree(os.path.dirname(workDir))
- os.makedirs(os.path.dirname(workDir),exist_ok=True)
-
-if __name__ == '__main__':
- main()
+#!/ usr / bin / env python3
+
+import os import shutil
+
+ def main() : nofJobs = 500 timeLimit = "08:00:00" geoSetup = "sis100_electron" plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
+
+#ADDED
+ mirrorRotation = 15 dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/" jobName = "RICH" + str(mirrorRotation)
+
+#All data in data dir will be removed
+ removeData = False
+#jobName = "LMVM"
+ if removeData : print("All data in dataDir will be removed. Dir:" + dataDir) print("Removing...") if os.path.exists(dataDir) : shutil.rmtree(dataDir) os.makedirs(dataDir, exist_ok = True)
+
+ for plutoParticle in plutoParticles : dataDirPluto = dataDir + plutoParticle logFile = dataDirPluto + "/log/log_slurm-%A_%a.out" errorFile = dataDirPluto + "/error/error_slurm-%A_%a.out" workDir = dataDirPluto + "/workdir/" makeLogErrorDirs(logFile, errorFile, workDir)
+
+#- p debug
+ commandStr =('sbatch --job-name={} --time={} --output={} --error={} --array=1-{} -- ./batch_job.py {} {} {}').format(jobName, timeLimit, logFile, errorFile, nofJobs, dataDirPluto, geoSetup, plutoParticle)
+
+ os.system(commandStr)
+
+
+ def makeLogErrorDirs(logFile, errorFile, workDir) : if os.path.exists(os.path.dirname(logFile)) : shutil.rmtree(os.path.dirname(logFile)) os.makedirs(os.path.dirname(logFile), exist_ok = True)
+
+ if os.path.exists(os.path.dirname(errorFile)) : shutil.rmtree(os.path.dirname(errorFile)) os.makedirs(os.path.dirname(errorFile), exist_ok = True)
+
+ if os.path.exists(os.path.dirname(workDir)) : shutil.rmtree(os.path.dirname(workDir)) os.makedirs(os.path.dirname(workDir), exist_ok = True)
+
+ if __name__ == '__main__' : main()
diff --git a/macro/analysis/dielectron/batch_send_common.py b/macro/analysis/dielectron/batch_send_common.py
new file mode 100644
index 0000000000000000000000000000000000000000..288932097ac2315195ed7f8cb4c267c178069c69
--- /dev/null
+++ b/macro/analysis/dielectron/batch_send_common.py
@@ -0,0 +1,56 @@
+#!/ usr / bin / env python3
+
+## #changed on Feb 8th; first check if works or not ###
+
+import os
+import shutil
+
+def main():
+ nofJobs = 100
+ timeLimit="08:00:00"
+ geoSetup="sis100_electron"
+ plutoParticles = ["inmed_had_epem", "wdalitz", "w", "phi", "qgp_epem"]
+
+ recoDataDirIn = "/lustre/cbm/pwg/common/mc/cbmsim/apr20_fr_18.2.1_fs_jun19p1/urqmd_pluto_"
+ dataDirOut = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm_CommonData/"
+
+ jobName = "LMVM"
+
+#All data in data dir will be removed
+ removeData = False
+ if removeData:
+ print("All data in dataDirOut will be removed. Dir:" + dataDirOut)
+ print("Removing...")
+ if os.path.exists(dataDirOut):
+ shutil.rmtree(dataDirOut)
+ os.makedirs(dataDirOut,exist_ok=True)
+
+ for plutoParticle in plutoParticles:
+ recoDataDirInPluto = recoDataDirIn + plutoParticle + "/"
+ dataDirOutPluto = dataDirOut + plutoParticle
+ logFile = dataDirOutPluto + "/log/log_slurm-%A_%a.out"
+ errorFile = dataDirOutPluto + "/error/error_slurm-%A_%a.out"
+ workDir = dataDirOutPluto + "/workdir/"
+ makeLogErrorDirs(logFile, errorFile, workDir)
+
+#- p debug
+ commandStr=('sbatch --job-name={} --time={} --output={} --error={} --array=1-{} -- ./batch_job_common.py {} {} {} {}').format(jobName, timeLimit, logFile, errorFile, nofJobs, recoDataDirInPluto, dataDirOutPluto, geoSetup, plutoParticle)
+
+ os.system(commandStr)
+
+
+def makeLogErrorDirs(logFile, errorFile, workDir):
+ if os.path.exists(os.path.dirname(logFile)):
+ shutil.rmtree(os.path.dirname(logFile))
+ os.makedirs(os.path.dirname(logFile),exist_ok=True)
+
+ if os.path.exists(os.path.dirname(errorFile)):
+ shutil.rmtree(os.path.dirname(errorFile))
+ os.makedirs(os.path.dirname(errorFile),exist_ok=True)
+
+ if os.path.exists(os.path.dirname(workDir)):
+ shutil.rmtree(os.path.dirname(workDir))
+ os.makedirs(os.path.dirname(workDir),exist_ok=True)
+
+if __name__ == '__main__':
+ main()
diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py
index bc0ea164b0ef3a47179b6af3449fd8ef66e16056..2de2035e4b349df2b192af31f6be8c431bf9dad2 100644
--- a/macro/analysis/dielectron/draw_all.py
+++ b/macro/analysis/dielectron/draw_all.py
@@ -1,32 +1,20 @@
-#!/usr/bin/env python3
+#!/ usr / bin / env python3
import os, shutil
-def main():
- plutoParticles = ["inmed", "phi", "omegaepem", "omegadalitz", "qgp"]
- dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/"
- dataDirHistos = dataDir + "histograms/"
- macroDir = "/lustre/cbm/users/criesen/cbmroot/macro/analysis/dielectron"
-
- if os.path.exists(dataDirHistos):
- shutil.rmtree(dataDirHistos)
- os.mkdir(dataDirHistos)
-
- for plutoParticle in plutoParticles:
-
- outFile = dataDirHistos + plutoParticle
-
- inFilesAna = dataDir + plutoParticle + "/analysis.all.root"
- uMvd = False
- drawSig = True
- os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesAna, outFile, uMvd, drawSig))
-
- inFilesLitqa = dataDir + plutoParticle + "/litqa.all.root"
- os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesLitqa, outFile))
-
-
- print("======================================================================================= DRAW ALL ==================================================================================")
- os.system(('root -l -b -q {}/draw_analysis_all.C').format(macroDir)) # default values of macro will be taken (see there)
-
-if __name__ == '__main__':
- main()
+ def main() : plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"] dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/" dataDirHistos = dataDir + "histograms/" macroDir = "/lustre/cbm/users/criesen/cbmroot/macro/analysis/dielectron"
+
+ if os.path.exists(dataDirHistos) : shutil.rmtree(dataDirHistos) os.mkdir(dataDirHistos)
+
+ for plutoParticle in plutoParticles :
+
+ outFile = dataDirHistos + plutoParticle
+
+ inFilesAna = dataDir + plutoParticle + "/analysis.all.root" uMvd = True drawSig = True os.system(('root -l -b -q {}/draw_analysis.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesAna, outFile, uMvd, drawSig))
+
+ inFilesLitqa = dataDir + plutoParticle + "/litqa.all.root" os.system(('root -l -b -q {}/draw_litqa.C\(\\"{}\\",\\"{}\\"\)').format(macroDir, inFilesLitqa, outFile))
+
+
+ print("===== DRAW ALL =====") os.system(('root -l -b -q {}/draw_analysis_all.C').format(macroDir))
+
+ if __name__ == '__main__' : main()
diff --git a/macro/analysis/dielectron/draw_litqa.C b/macro/analysis/dielectron/draw_litqa.C
index dd730e935b0832631d30ffb4d985025282e224bc..df0eaaa77f26787b751fb4174558e229ca695f9d 100644
--- a/macro/analysis/dielectron/draw_litqa.C
+++ b/macro/analysis/dielectron/draw_litqa.C
@@ -12,11 +12,11 @@ void draw_litqa(const string& histRootFile = "/lustre/nyx/cbm/users/criesen/cbm/
CbmSimulationReport* trackingQaReport = new CbmLitTrackingQaReport();
trackingQaReport->Create(histRootFile, outputDirTracking);
- // string outputDirClustering = resultDir + "results_litqa_clustering/";
- // gSystem->mkdir(outputDirClustering.c_str(), true);
+ string outputDirClustering = resultDir + "results_litqa_clustering/";
+ gSystem->mkdir(outputDirClustering.c_str(), true);
- // CbmSimulationReport* clusteringQaReport = new CbmLitClusteringQaReport();
- // clusteringQaReport->Create(fileName, outputDirClustering);
+ CbmSimulationReport* clusteringQaReport = new CbmLitClusteringQaReport();
+ clusteringQaReport->Create(histRootFile, outputDirClustering);
// CbmSimulationReport* fitQaReport = new CbmLitFitQaReport();
// fitQaReport->Create(fileName, outputDir);
diff --git a/macro/analysis/dielectron/hadd_many.py b/macro/analysis/dielectron/hadd_many.py
old mode 100755
new mode 100644
index f22dc824020c0dac504c6d6b1cfcd958786c6a9e..249f3357d1f4045141d9e76559942aca933bc70f
--- a/macro/analysis/dielectron/hadd_many.py
+++ b/macro/analysis/dielectron/hadd_many.py
@@ -1,24 +1,11 @@
-import os
-import shutil
+import os import shutil
-def main():
- plutoParticles = ["inmed", "phi", "omegaepem", "omegadalitz", "qgp"];
- dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm"
-
- for plutoParticle in plutoParticles:
- dataDirPluto = dataDir + "/" + plutoParticle
+ def main() : plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"] dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm"
- inFilesLitQa = dataDirPluto + "/litqa.*.root"
- outFileLitQa = dataDirPluto + "/litqa.all.root"
- if os.path.exists(outFileLitQa):
- os.remove(outFileLitQa)
- os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa))
+ for plutoParticle in plutoParticles : dataDirPluto = dataDir + "/" + plutoParticle
- inFilesAna = dataDirPluto + "/analysis.*.root"
- outFileAna = dataDirPluto + "/analysis.all.root"
- if os.path.exists(outFileAna):
- os.remove(outFileAna)
- os.system(('hadd -j -T -f {} {}').format(outFileAna, inFilesAna))
+ inFilesLitQa = dataDirPluto + "/litqa.*.root" outFileLitQa = dataDirPluto + "/litqa.all.root" if os.path.exists(outFileLitQa) : os.remove(outFileLitQa) os.system(('hadd -j -T -f {} {}').format(outFileLitQa, inFilesLitQa))
-if __name__ == '__main__':
- main()
+ inFilesAna = dataDirPluto + "/analysis.*.root" outFileAna = dataDirPluto + "/analysis.all.root" if os.path.exists(outFileAna) : os.remove(outFileAna) os.system(('hadd -j -T -f {} {}').format(outFileAna, inFilesAna))
+
+ if __name__ == '__main__' : main()
diff --git a/macro/analysis/dielectron/run_litqa.C b/macro/analysis/dielectron/run_litqa.C
index 44f61e53055ae71c2161ad014f7f0be05a15e287..8808346e0f351c2685da3515ac3a0d597ffb1364 100644
--- a/macro/analysis/dielectron/run_litqa.C
+++ b/macro/analysis/dielectron/run_litqa.C
@@ -34,9 +34,9 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/
// - TOF digitisation parameters
if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
- TObjString* tofFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
- parFileList->Add(tofFile);
- std::cout << "-I- " << myName << ": Using parameter file " << tofFile->GetString() << std::endl;
+ //TObjString* tofFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
+ //parFileList->Add(tofFile);
+ //std::cout << "-I- " << myName << ": Using parameter file " << tofFile->GetString() << std::endl;
TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
parFileList->Add(tofBdfFile);
std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
@@ -93,7 +93,7 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/
run->AddTask(trackingQa);
CbmLitFitQa* fitQa = new CbmLitFitQa();
- fitQa->SetMvdMinNofHits(0);
+ fitQa->SetMvdMinNofHits(3);
fitQa->SetStsMinNofHits(6);
fitQa->SetMuchMinNofHits(10);
fitQa->SetTrdMinNofHits(2);
@@ -102,7 +102,7 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/
CbmLitClusteringQa* clusteringQa = new CbmLitClusteringQa();
clusteringQa->SetOutputDir("");
- // run->AddTask(clusteringQa);
+ run->AddTask(clusteringQa);
CbmLitTofQa* tofQa = new CbmLitTofQa();
tofQa->SetOutputDir(std::string(""));
diff --git a/macro/analysis/dielectron/run_reco.C b/macro/analysis/dielectron/run_reco.C
index 5dcd2a24da896b0fc9922da487cd758e24188e9e..796fd207c69c095883b38a54adbdf2751a07074f 100644
--- a/macro/analysis/dielectron/run_reco.C
+++ b/macro/analysis/dielectron/run_reco.C
@@ -34,9 +34,9 @@ void run_reco(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/l
// - TOF digitisation parameters
if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
- TObjString* tofFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
- parFileList->Add(tofFile);
- std::cout << "-I- " << myName << ": Using parameter file " << tofFile->GetString() << std::endl;
+ //TObjString* tofFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par");
+ //parFileList->Add(tofFile);
+ //std::cout << "-I- " << myName << ": Using parameter file " << tofFile->GetString() << std::endl;
TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
parFileList->Add(tofBdfFile);
std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
diff --git a/reco/littrack/cbm/qa/tracking/CbmLitTrackingQaReport.cxx b/reco/littrack/cbm/qa/tracking/CbmLitTrackingQaReport.cxx
index 6343c420d0a4d65d1d3b6d2c704203e07696fdfb..be0191445fcc03bd9fd645fcef5eef8068f656bc 100644
--- a/reco/littrack/cbm/qa/tracking/CbmLitTrackingQaReport.cxx
+++ b/reco/littrack/cbm/qa/tracking/CbmLitTrackingQaReport.cxx
@@ -541,8 +541,8 @@ void CbmLitTrackingQaReport::CalculatePionSuppressionHistos()
string recHistName = FindAndReplace(psHistName, "_PionSup_", "_Rec_");
string pionRecHistName = FindAndReplace(psHistName, "_PionSup_", "_RecPions_");
DivideHistos(HM()->H1(pionRecHistName), HM()->H1(recHistName), psHist, 1.);
- // psHist->SetMinimum(1.);
- // psHist->SetMaximum(20000.);
+ psHist->SetMinimum(10.);
+ psHist->SetMaximum(2e5);
}
}