From cc53d8e1d83797ad6fd80772554daf8f47e7b363 Mon Sep 17 00:00:00 2001
From: drslebedev <dr.s.lebedev@gmail.com>
Date: Wed, 30 Jun 2021 12:41:47 +0200
Subject: [PATCH] Update lmvm python scripts
Adding lmvm combinatorial pairs procedure
Develop combinatorial BG procedure and add corresponding histograms
Correct pairing in CB and reverse formatting of python macros
Add Single Particle Yield Histos
Particle Yield and Acceptance of Electrons and Positrons is observed for various detector combinations.
Formatting of python scripts after they had been tored by wrong application of clang format.
WIP Add single particle histos
WIP: Update single particle histograms
WIP: Add histograms for nPoint
Number of points the electrons / positrons leave in the various detectors
Modify macros for common production
Also some small adaptions in CbmAnaDielectronTask.cxx and in ana and litqa macros were done.
Apply clang format
---
.../dielectron/lmvm/CbmAnaDielectronTask.cxx | 348 ++++++++++++------
.../dielectron/lmvm/CbmAnaDielectronTask.h | 22 +-
.../lmvm/CbmAnaDielectronTaskDraw.cxx | 325 ++++++++++++++--
.../lmvm/CbmAnaDielectronTaskDrawAll.cxx | 256 ++++++-------
macro/analysis/dielectron/batch_job.py | 86 +++--
macro/analysis/dielectron/batch_job_common.py | 146 +++++---
macro/analysis/dielectron/batch_send.py | 82 +++--
.../analysis/dielectron/batch_send_common.py | 12 +-
macro/analysis/dielectron/draw_all.py | 43 ++-
macro/analysis/dielectron/draw_analysis.C | 0
macro/analysis/dielectron/draw_analysis_all.C | 0
macro/analysis/dielectron/hadd_many.py | 33 +-
macro/analysis/dielectron/run_analysis.C | 18 +-
macro/analysis/dielectron/run_litqa.C | 2 +-
14 files changed, 942 insertions(+), 431 deletions(-)
mode change 100644 => 100755 macro/analysis/dielectron/batch_job.py
mode change 100644 => 100755 macro/analysis/dielectron/batch_job_common.py
mode change 100644 => 100755 macro/analysis/dielectron/batch_send.py
mode change 100644 => 100755 macro/analysis/dielectron/batch_send_common.py
mode change 100644 => 100755 macro/analysis/dielectron/draw_all.py
mode change 100644 => 100755 macro/analysis/dielectron/draw_analysis.C
mode change 100644 => 100755 macro/analysis/dielectron/draw_analysis_all.C
mode change 100644 => 100755 macro/analysis/dielectron/hadd_many.py
mode change 100644 => 100755 macro/analysis/dielectron/run_analysis.C
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
index e5fb3add05..677a3686e2 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx
@@ -69,85 +69,6 @@
#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);
@@ -279,6 +200,12 @@ CbmAnaDielectronTask::CbmAnaDielectronTask()
, fh_nof_plutoPositrons()
, fh_nof_urqmdElectrons()
, fh_nof_urqmdPositrons()
+ , fh_nof_plutoElectrons_p_pt()
+ , fh_nof_plutoPositrons_p_pt()
+ , fh_nof_urqmdElectrons_p_pt()
+ , fh_nof_urqmdPositrons_p_pt()
+ , fh_nof_particles_acc()
+ , fh_nof_points()
, fh_pi0_minv()
, fh_eta_minv()
, fh_gamma_minv()
@@ -547,14 +474,24 @@ void CbmAnaDielectronTask::InitHists()
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,
+ // 1D Histograms for single particle yield vs. momentum
+ CreateAnalysisStepsH1(fh_nof_plutoElectrons, "fh_nof_plutoElectrons", "P [Gev/c]", "yield pluto electrons", 100, 0,
10.);
- CreateAnalysisStepsH1(fh_nof_plutoPositrons, "fh_nof_plutoPositrons", "P [Gev/c]", "number pluto positrons", 10000, 0,
+ CreateAnalysisStepsH1(fh_nof_plutoPositrons, "fh_nof_plutoPositrons", "P [Gev/c]", "yield pluto positrons", 100, 0,
10.);
- CreateAnalysisStepsH1(fh_nof_urqmdElectrons, "fh_nof_urqmdElectrons", "P [Gev/c]", "number urqmd electrons", 10000, 0,
+ CreateAnalysisStepsH1(fh_nof_urqmdElectrons, "fh_nof_urqmdElectrons", "P [Gev/c]", "yield urqmd electrons", 100, 0,
10.);
- CreateAnalysisStepsH1(fh_nof_urqmdPositrons, "fh_nof_urqmdPositrons", "P [Gev/c]", "number urqmd positrons", 10000, 0,
+ CreateAnalysisStepsH1(fh_nof_urqmdPositrons, "fh_nof_urqmdPositrons", "P [Gev/c]", "yield urqmd positrons", 100, 0,
10.);
+ // 2D Histograms for single particle yield vs. momentum and Pt
+ CreateAnalysisStepsH2(fh_nof_plutoElectrons_p_pt, "fh_nof_plutoElectrons_p_pt", "P [GeV/c]", "P_{t} [Gev/c]",
+ "yield pluto electrons", 100, 0, 10., 40, 0., 4.);
+ CreateAnalysisStepsH2(fh_nof_plutoPositrons_p_pt, "fh_nof_plutoPositrons_p_pt", "P [GeV/c]", "P_{t} [Gev/c]",
+ "yield pluto positrons", 100, 0, 10., 40, 0., 4.);
+ CreateAnalysisStepsH2(fh_nof_urqmdElectrons_p_pt, "fh_nof_urqmdElectrons_p_pt", "P [GeV/c]", "P_{t} [Gev/c]",
+ "yield urqmd electrons", 100, 0, 10., 40, 0., 4.);
+ CreateAnalysisStepsH2(fh_nof_urqmdPositrons_p_pt, "fh_nof_urqmdPositrons_p_pt", "P [GeV/c]", "P_{t} [Gev/c]",
+ "yield urqmd positrons", 100, 0, 10., 40, 0., 4.);
// 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.);
@@ -685,6 +622,63 @@ void CbmAnaDielectronTask::InitHists()
fHistoList.push_back(fh_mom_likelihood_El);
fHistoList.push_back(fh_mom_likelihood_Pi);
+
+ // Acceptance of single particles vs. momentum for various detector combinations
+ fh_nof_particles_acc.resize(20);
+ fh_nof_particles_acc[0] = new TH1D("fh_nof_particles_acc_pEl_mc", "fh_nof_particles_acc_pEl_mc", 100., 0., 10.);
+ fh_nof_particles_acc[1] = new TH1D("fh_nof_particles_acc_pPos_mc", "fh_nof_particles_acc_pPos_mc", 100., 0., 10.);
+ fh_nof_particles_acc[2] = new TH1D("fh_nof_particles_acc_pEl_sts", "fh_nof_particles_acc_pEl_sts", 100., 0., 10.);
+ fh_nof_particles_acc[3] = new TH1D("fh_nof_particles_acc_pPos_sts", "fh_nof_particles_acc_pPos_sts", 100., 0., 10.);
+ fh_nof_particles_acc[4] = new TH1D("fh_nof_particles_acc_pEl_rich", "fh_nof_particles_acc_pEl_rich", 100., 0., 10.);
+ fh_nof_particles_acc[5] = new TH1D("fh_nof_particles_acc_pPos_rich", "fh_nof_particles_acc_pPos_rich", 100., 0., 10.);
+ fh_nof_particles_acc[6] = new TH1D("fh_nof_particles_acc_pEl_trd", "fh_nof_particles_acc_pEl_trd", 100., 0., 10.);
+ fh_nof_particles_acc[7] = new TH1D("fh_nof_particles_acc_pPos_trd", "fh_nof_particles_acc_pPos_trd", 100., 0., 10.);
+ fh_nof_particles_acc[8] = new TH1D("fh_nof_particles_acc_pEl_tof", "fh_nof_particles_acc_pEl_tof", 100., 0., 10.);
+ fh_nof_particles_acc[9] = new TH1D("fh_nof_particles_acc_pPos_tof", "fh_nof_particles_acc_pPos_tof", 100., 0., 10.);
+ fh_nof_particles_acc[10] = new TH1D("fh_nof_particles_acc_uEl_mc", "fh_nof_particles_acc_uEl_mc", 100., 0., 10.);
+ fh_nof_particles_acc[11] = new TH1D("fh_nof_particles_acc_uPos_mc", "fh_nof_particles_acc_uPos_mc", 100., 0., 10.);
+ fh_nof_particles_acc[12] = new TH1D("fh_nof_particles_acc_uEl_sts", "fh_nof_particles_acc_uEl_sts", 100., 0., 10.);
+ fh_nof_particles_acc[13] = new TH1D("fh_nof_particles_acc_uPos_sts", "fh_nof_particles_acc_uPos_sts", 100., 0., 10.);
+ fh_nof_particles_acc[14] = new TH1D("fh_nof_particles_acc_uEl_rich", "fh_nof_particles_acc_uEl_rich", 100., 0., 10.);
+ fh_nof_particles_acc[15] =
+ new TH1D("fh_nof_particles_acc_uPos_rich", "fh_nof_particles_acc_uPos_rich", 100., 0., 10.);
+ fh_nof_particles_acc[16] = new TH1D("fh_nof_particles_acc_uEl_trd", "fh_nof_particles_acc_uEl_trd", 100., 0., 10.);
+ fh_nof_particles_acc[17] = new TH1D("fh_nof_particles_acc_uPos_trd", "fh_nof_particles_acc_uPos_trd", 100., 0., 10.);
+ fh_nof_particles_acc[18] = new TH1D("fh_nof_particles_acc_uEl_tof", "fh_nof_particles_acc_uEl_tof", 100., 0., 10.);
+ fh_nof_particles_acc[19] = new TH1D("fh_nof_particles_acc_uPos_tof", "fh_nof_particles_acc_uPos_tof", 100., 0., 10.);
+
+ int nHist = fh_nof_particles_acc.size();
+ for (Int_t i = 0; i < nHist; i++) {
+ fh_nof_particles_acc[i]->GetXaxis()->SetTitle("P [GeV/c]");
+ fh_nof_particles_acc[i]->GetYaxis()->SetTitle("Yield");
+ fHistoList.push_back(fh_nof_particles_acc[i]);
+ }
+
+ // Number of points the electrons and positrons left in various detectors
+ fh_nof_points.resize(16);
+ fh_nof_points[0] = new TH1D("fh_nof_points_pEl_sts", "fh_nof_points_pEl_sts", 50., 0., 50.);
+ fh_nof_points[1] = new TH1D("fh_nof_points_pPos_sts", "fh_nof_points_pPos_sts", 50., 0., 50.);
+ fh_nof_points[2] = new TH1D("fh_nof_points_pEl_rich", "fh_nof_points_pEl_rich", 50., 0., 50.);
+ fh_nof_points[3] = new TH1D("fh_nof_points_pPos_rich", "fh_nof_points_pPos_rich", 50., 0., 50.);
+ fh_nof_points[4] = new TH1D("fh_nof_points_pEl_trd", "fh_nof_points_pEl_trd", 50., 0., 50.);
+ fh_nof_points[5] = new TH1D("fh_nof_points_pPos_trd", "fh_nof_points_pPos_trd", 50., 0., 50.);
+ fh_nof_points[6] = new TH1D("fh_nof_points_pEl_tof", "fh_nof_points_pEl_tof", 50., 0., 50.);
+ fh_nof_points[7] = new TH1D("fh_nof_points_pPos_tof", "fh_nof_points_pPos_tof", 50., 0., 50.);
+ fh_nof_points[8] = new TH1D("fh_nof_points_uEl_sts", "fh_nof_points_uEl_sts", 50., 0., 50.);
+ fh_nof_points[9] = new TH1D("fh_nof_points_uPos_sts", "fh_nof_points_uPos_sts", 50., 0., 50.);
+ fh_nof_points[10] = new TH1D("fh_nof_points_uEl_rich", "fh_nof_points_uEl_rich", 50., 0., 50.);
+ fh_nof_points[11] = new TH1D("fh_nof_points_uPos_rich", "fh_nof_points_uPos_rich", 50., 0., 50.);
+ fh_nof_points[12] = new TH1D("fh_nof_points_uEl_trd", "fh_nof_points_uEl_trd", 50., 0., 50.);
+ fh_nof_points[13] = new TH1D("fh_nof_points_uPos_trd", "fh_nof_points_uPos_trd", 50., 0., 50.);
+ fh_nof_points[14] = new TH1D("fh_nof_points_uEl_tof", "fh_nof_points_uEl_tof", 50., 0., 50.);
+ fh_nof_points[15] = new TH1D("fh_nof_points_uPos_tof", "fh_nof_points_uPos_tof", 50., 0., 50.);
+
+ int nHistPoints = fh_nof_points.size();
+ for (Int_t i = 0; i < nHistPoints; i++) {
+ fh_nof_points[i]->GetXaxis()->SetTitle("nofPoints");
+ fh_nof_points[i]->GetYaxis()->SetTitle("Entries");
+ fHistoList.push_back(fh_nof_points[i]);
+ }
}
InitStatus CbmAnaDielectronTask::Init()
@@ -805,9 +799,8 @@ void CbmAnaDielectronTask::Exec(Option_t*)
else {
Fatal("CbmAnaDielectronTask::Exec", "No PrimaryVertex array!");
}
- // cout << "I'm here" << endl;
// CbmLitMCTrackCreator::Instance()->CreateReco();
- //
+
if (useMbias || (!useMbias && isCentralCollision)) {
FillRichRingNofHits();
MCPairs();
@@ -985,52 +978,153 @@ void CbmAnaDielectronTask::FillNofChargedParticles()
cout << "FillNofChargedParticles: nofMcTracks = " << nofMcTracks << endl;
cout << "FillNofChargedParticles: nCand = " << nCand << endl;
+ int nChargeMC = 0;
+ int nChargeCand = 0;
+
for (Int_t i = 0; i < nofMcTracks; i++) {
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 momT = mcTrack->GetPt();
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 isMcTrAcc =
- IsMcTrackAccepted(i); // before: = (mcTrack->GetNPoints(ECbmModuleId::kSts) >= 4) ? true : false;
+ Bool_t isUrqmdElectron = (!isMcSignalElectron && isMcElectron) ? true : false;
+ Bool_t isMcTrAcc = IsMcTrackAccepted(i);
Bool_t isPrimary = (motherId == -1) ? true : false;
- if (charge != 0) cout << "charge = " << charge << endl;
+ CbmMCTrack* tr = (CbmMCTrack*) fMCTracks->At(i);
+ Int_t nRichPoints = fNofHitsInRingMap[i];
+ int nStsPoints = tr->GetNPoints(ECbmModuleId::kSts);
+ int nTrdPoints = tr->GetNPoints(ECbmModuleId::kTrd);
+ int nTofPoints = tr->GetNPoints(ECbmModuleId::kTof);
+
+ bool isStsAcc = (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4) ? true : false;
+ bool isRichAcc =
+ (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4 && nRichPoints >= 7) ? true : false;
+ bool isTrdAcc = (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4 && nRichPoints >= 7
+ && tr->GetNPoints(ECbmModuleId::kTrd) >= 2)
+ ? true
+ : false;
+ bool isTofAcc = (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4 && nRichPoints >= 7
+ && tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 1)
+ ? true
+ : false; // Mind: kTof was '0' before! changed on 2.7.21; same in IsMcTrackAccepted(int)
+
+ if (charge != 0) {
+ cout << "charge (from GetCharge()) = " << charge << endl;
+ nChargeMC++;
+ }
if (!isMcSignalElectron && isMcTrackCharged && isPrimary) nofChargedUrqmdParticles++;
if (!isMcSignalElectron && isMcTrackCharged && isMcTrAcc && isPrimary) nofChargedUrqmdParticlesAcc++;
+ // 1D Histos: Yield vs. Momentum
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);
+ if (isUrqmdElectron && charge < 0) fh_nof_urqmdElectrons[kMc]->Fill(mom);
+ if (isUrqmdElectron && charge > 0) fh_nof_urqmdPositrons[kMc]->Fill(mom);
+ if (isUrqmdElectron && isMcTrAcc && charge < 0) fh_nof_urqmdElectrons[kAcc]->Fill(mom);
+ if (isUrqmdElectron && isMcTrAcc && charge > 0) fh_nof_urqmdPositrons[kAcc]->Fill(mom);
+
+ // 2D Histos: Yield vs. Momentum and Pt
+ if (isMcSignalElectron && charge < 0) fh_nof_plutoElectrons_p_pt[kMc]->Fill(mom, momT);
+ if (isMcSignalElectron && charge > 0) fh_nof_plutoPositrons_p_pt[kMc]->Fill(mom, momT);
+ if (isMcSignalElectron && isMcTrAcc && charge < 0) fh_nof_plutoElectrons_p_pt[kAcc]->Fill(mom, momT);
+ if (isMcSignalElectron && isMcTrAcc && charge > 0) fh_nof_plutoPositrons_p_pt[kAcc]->Fill(mom, momT);
+ if (isUrqmdElectron && charge < 0) fh_nof_urqmdElectrons_p_pt[kMc]->Fill(mom, momT);
+ if (isUrqmdElectron && charge > 0) fh_nof_urqmdPositrons_p_pt[kMc]->Fill(mom, momT);
+ if (isUrqmdElectron && isMcTrAcc && charge < 0) fh_nof_urqmdElectrons_p_pt[kAcc]->Fill(mom, momT);
+ if (isUrqmdElectron && isMcTrAcc && charge > 0) fh_nof_urqmdPositrons_p_pt[kAcc]->Fill(mom, momT);
+
+ // Checking Acceptance of diff. Detectors
+ if (isMcSignalElectron && charge < 0) fh_nof_particles_acc[0]->Fill(mom); // PLUTO electron
+ if (isMcSignalElectron && charge > 0) fh_nof_particles_acc[1]->Fill(mom); // PLUTO positron
+ if (isMcSignalElectron && isStsAcc && charge < 0) fh_nof_particles_acc[2]->Fill(mom); // PLUTO electron
+ if (isMcSignalElectron && isStsAcc && charge > 0) fh_nof_particles_acc[3]->Fill(mom); // PLUTO positron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && charge < 0)
+ fh_nof_particles_acc[4]->Fill(mom); // PLUTO electron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && charge > 0)
+ fh_nof_particles_acc[5]->Fill(mom); // PLUTO positron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && isTrdAcc && charge < 0)
+ fh_nof_particles_acc[6]->Fill(mom); // PLUTO electron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && isTrdAcc && charge > 0)
+ fh_nof_particles_acc[7]->Fill(mom); // PLUTO positron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && isTrdAcc && isTofAcc && charge < 0)
+ fh_nof_particles_acc[8]->Fill(mom); // PLUTO electron
+ if (isMcSignalElectron && isStsAcc && isRichAcc && isTrdAcc && isTofAcc && charge > 0)
+ fh_nof_particles_acc[9]->Fill(mom); // PLUTO positron
+ if (isUrqmdElectron && charge < 0) fh_nof_particles_acc[10]->Fill(mom); // UrQMD electron
+ if (isUrqmdElectron && charge > 0) fh_nof_particles_acc[11]->Fill(mom); // UrQMD positron
+ if (isUrqmdElectron && isStsAcc && charge < 0) fh_nof_particles_acc[12]->Fill(mom); // UrQMD electron
+ if (isUrqmdElectron && isStsAcc && charge > 0) fh_nof_particles_acc[13]->Fill(mom); // UrQMD positron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && charge < 0) fh_nof_particles_acc[14]->Fill(mom); // UrQMD electron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && charge > 0) fh_nof_particles_acc[15]->Fill(mom); // UrQMD positron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && isTrdAcc && charge < 0)
+ fh_nof_particles_acc[16]->Fill(mom); // UrQMD electron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && isTrdAcc && charge > 0)
+ fh_nof_particles_acc[17]->Fill(mom); // UrQMD positron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && isTrdAcc && isTofAcc && charge < 0)
+ fh_nof_particles_acc[18]->Fill(mom); // UrQMD electron
+ if (isUrqmdElectron && isStsAcc && isRichAcc && isTrdAcc && isTofAcc && charge > 0)
+ fh_nof_particles_acc[19]->Fill(mom); // UrQMD positron
+
+ // Fill histos with number of points
+ if (isMcSignalElectron && charge < 0) {
+ fh_nof_points[0]->Fill(nStsPoints);
+ fh_nof_points[2]->Fill(nRichPoints);
+ fh_nof_points[4]->Fill(nTrdPoints);
+ fh_nof_points[6]->Fill(nTofPoints);
+ }
+ if (isMcSignalElectron && charge > 0) {
+ fh_nof_points[1]->Fill(nStsPoints);
+ fh_nof_points[3]->Fill(nRichPoints);
+ fh_nof_points[5]->Fill(nTrdPoints);
+ fh_nof_points[7]->Fill(nTofPoints);
+ }
+ if (isUrqmdElectron && charge < 0) {
+ fh_nof_points[8]->Fill(nStsPoints);
+ fh_nof_points[10]->Fill(nRichPoints);
+ fh_nof_points[12]->Fill(nTrdPoints);
+ fh_nof_points[14]->Fill(nTofPoints);
+ }
+ if (isUrqmdElectron && charge > 0) {
+ fh_nof_points[9]->Fill(nStsPoints);
+ fh_nof_points[11]->Fill(nRichPoints);
+ fh_nof_points[13]->Fill(nTrdPoints);
+ fh_nof_points[15]->Fill(nTofPoints);
+ }
}
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;
+ int charge = fCandidates[i].fCharge;
+ if (charge != 0) {
+ cout << "fCandidates[i].fCharge = " << charge << endl;
+ nChargeCand++;
+ }
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;
+ Bool_t isMcElectron =
+ (fCandidates[i].fIsMcPi0Electron || fCandidates[i].fIsMcEtaElectron || fCandidates[i].fIsMcGammaElectron) ? true
+ : false;
+ Bool_t isUrqmdElectron = (!isMcSignalElectron && isMcElectron) ? true : false;
+ Bool_t isChiPrimary = fCandidates[i].fChi2Prim < fCuts.fChiPrimCut;
+ Bool_t isElectron = fCandidates[i].fIsElectron;
+ Bool_t isGammaCut = !fCandidates[i].fIsGamma;
+ //Bool_t isMvd1Cut = fCandidates[i].fIsMvd1CutElectron;
+ //Bool_t isMvd2Cut = fCandidates[i].fIsMvd2CutElectron;
+ Bool_t isStCut = fCandidates[i].fIsStCutElectron;
+ Bool_t isRtCut = fCandidates[i].fIsRtCutElectron;
+ Bool_t isTtCut = fCandidates[i].fIsTtCutElectron;
+ Bool_t isPtCut = fCandidates[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());
@@ -1038,6 +1132,13 @@ void CbmAnaDielectronTask::FillNofChargedParticles()
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 (isChiPrimary && isElectron && charge < 0)
+ fh_nof_plutoElectrons[kChi2Prim]->Fill(
+ fCandidates[i]
+ .fMomentum.Mag()); // !!!TODO: is NOT kChi2Prim step but El-ID wo. MC info; to show difference between with
+ if (isChiPrimary && isElectron && charge > 0)
+ fh_nof_plutoPositrons[kChi2Prim]->Fill(
+ fCandidates[i].fMomentum.Mag()); // ... and wo MC info!! Remove afterwards !!
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)
@@ -1049,23 +1150,31 @@ void CbmAnaDielectronTask::FillNofChargedParticles()
&& 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)
+ if (isUrqmdElectron && charge < 0) fh_nof_urqmdElectrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isUrqmdElectron && charge > 0) fh_nof_urqmdPositrons[kReco]->Fill(fCandidates[i].fMomentum.Mag());
+ if (isUrqmdElectron && isChiPrimary && isElectron && charge < 0)
fh_nof_urqmdElectrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
- if (!isMcSignalElectron && isChiPrimary && isElectron && charge > 0)
+ if (isUrqmdElectron && isChiPrimary && isElectron && charge > 0)
fh_nof_urqmdPositrons[kElId]->Fill(fCandidates[i].fMomentum.Mag());
- if (!isMcSignalElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && charge < 0)
+ if (isChiPrimary && isElectron && charge < 0)
+ fh_nof_urqmdElectrons[kChi2Prim]->Fill(
+ fCandidates[i]
+ .fMomentum.Mag()); // !!!TODO: is NOT kChi2Prim step but El-ID wo. MC info; to show difference between with
+ if (isChiPrimary && isElectron && charge > 0)
+ fh_nof_urqmdPositrons[kChi2Prim]->Fill(
+ fCandidates[i].fMomentum.Mag()); // ... and wo MC info!! Remove afterwards !!
+ if (isUrqmdElectron && 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)
+ if (isUrqmdElectron && 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
+ if (isUrqmdElectron && 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
+ if (isUrqmdElectron && isChiPrimary && isElectron && isGammaCut && isStCut && isRtCut && isTtCut && isPtCut
&& charge > 0)
fh_nof_urqmdPositrons[kPtCut]->Fill(fCandidates[i].fMomentum.Mag());
}
+ cout << "nChargeMC = " << nChargeMC << ", nChargeCand = " << nChargeCand << endl;
}
Bool_t CbmAnaDielectronTask::IsMcTrackAccepted(Int_t mcTrackInd)
@@ -1074,7 +1183,7 @@ Bool_t CbmAnaDielectronTask::IsMcTrackAccepted(Int_t mcTrackInd)
if (tr == NULL) return false;
Int_t nRichPoints = fNofHitsInRingMap[mcTrackInd];
return (tr->GetNPoints(ECbmModuleId::kMvd) + tr->GetNPoints(ECbmModuleId::kSts) >= 4 && nRichPoints >= 7
- && tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 0);
+ && tr->GetNPoints(ECbmModuleId::kTrd) >= 2 && tr->GetNPoints(ECbmModuleId::kTof) > 1);
}
void CbmAnaDielectronTask::SingleParticleAcceptance()
@@ -1424,8 +1533,7 @@ void CbmAnaDielectronTask::CombinatorialPairs()
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;
+ double weight = (isSignal1 && isSignal2) ? fWeight : 1.;
CbmLmvmKinematicParams pRec =
CbmLmvmKinematicParams::KinematicParamsWithCandidates(&fCandidatesTotal[iCand1], &fCandidatesTotal[iCand2]);
@@ -2557,6 +2665,14 @@ void CbmAnaDielectronTask::Finish()
void CbmAnaDielectronTask::SetEnergyAndPlutoParticle(const string& energy, const string& particle)
{
+ // names of particles in common production differ from our names
+ if (particle == "rho0") particle == "inmed"; // or "inmed_had_epem"??
+ else if (particle == "wdalitz")
+ particle == "omegadalitz";
+ else if (particle == "w")
+ particle == "omegaepem";
+ else if (particle == "qgp_epem")
+ particle == "qgp";
// Au+Au centr old scaling factors
/* if (energy == "8gev" || energy == "10gev") {
@@ -2575,7 +2691,9 @@ void CbmAnaDielectronTask::SetEnergyAndPlutoParticle(const string& energy, const
//either old or new!!!
// Au+Au centr new scaling factors
}*/
- if (energy == "8gev") {
+ if (energy == "8gev"
+ || energy
+ == "12gev") { // TODO: 12 GeV was added only to check common production; must add weights for this energy
// weight omega = Multiplicity * Branching Ratio = 19 * 7.28e-5 for 8 AGeV beam energy
if (particle == "omegaepem") this->SetWeight(2.5 * 7.28e-5);
// weight omega = Multiplicity * Branching Ratio = 19 * 7.7e-4 for 8 AGeV beam energy
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
index d46ec889f4..7926e15075 100755
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.h
@@ -323,10 +323,23 @@ private:
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_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<TH2D*> fh_nof_plutoElectrons_p_pt; // nof electrons / positrons vs. momentum and Pt
+ std::vector<TH2D*> fh_nof_plutoPositrons_p_pt; // nof electrons / positrons vs. momentum and Pt
+ std::vector<TH2D*> fh_nof_urqmdElectrons_p_pt; // nof electrons / positrons vs. momentum and Pt
+ std::vector<TH2D*> fh_nof_urqmdPositrons_p_pt; // nof electrons / positrons vs. momentum and Pt
+
+ std::vector<TH1D*>
+ fh_nof_particles_acc; // To compare accepted with created PLUTO and UrQMD particles for various detector combinations
+
+ // Number of points the electrons and positrons left in various detectors
+ //[0]=Pluto-El STS, [1]=Pluto-Pos STS, [2]=Pluto-El RICH, [3]=Pluto-Pos RICH, [4]=Pluto-El TRD, [5]=Pluto-Pos TRD, [6]=Pluto-El ToF, [7]=Pluto-Pos ToF
+ //[8]=UrQMD-El STS, [9]=UrQMD-Pos STS, [10]=UrQMD-El RICH, [11]=UrQMD-Pos RICH, [12]=UrQMD-El TRD, [13]=UrQMD-Pos TRD, [14]=UrQMD-El ToF, [15]=UrQMD-Pos ToF
+ std::vector<TH1D*> fh_nof_points;
std::vector<TH1D*> fh_pi0_minv; // Invariant mass for Pi0
std::vector<TH1D*> fh_eta_minv; // Invariant mass for Eta
@@ -337,7 +350,6 @@ private:
std::vector<TH2D*> fh_eta_minv_pt; // Invariant mass vs. MC Pt
std::vector<TH2D*> fh_pi0_minv_pt; // Invariant mass vs. MC Pt
-
std::vector<TH1D*> fh_bg_truematch_minv; // Invariant mass for truly matched tracks
std::vector<TH1D*> fh_bg_truematch_el_minv; // Invariant mass for truly matched electron tracks
std::vector<TH1D*> fh_bg_truematch_notel_minv; // Invariant mass for truly matched tracks, not 2 electrons
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
index c3d3b24ad4..119013da61 100755
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDraw.cxx
@@ -120,6 +120,7 @@ void CbmAnaDielectronTaskDraw::DrawMomLikeHist()
void CbmAnaDielectronTaskDraw::DrawSingleParticleYield()
{
+ // Draw 1D Histos: Yield vs. Momentum
TCanvas* c2 = fHM->CreateCanvas("fh_mom_ElPos_pluto", "fh_mom_ElPos_pluto", 1000, 500);
c2->Divide(2, 1);
c2->cd(1);
@@ -238,24 +239,306 @@ void CbmAnaDielectronTaskDraw::DrawSingleParticleYield()
legendNofUP->AddEntry(nUrqmdElReco, "electrons kReco");
legendNofUP->AddEntry(nUrqmdPosReco, "positrons kReco");
legendNofUP->Draw();
+
+ // Acc. of single particle yield vs. momentum for various detector combinations
+ TH1D* nPlutoElMc2 = (TH1D*) H1("fh_nof_particles_acc_pEl_mc")->Clone();
+ TH1D* nPlutoPosMc2 = (TH1D*) H1("fh_nof_particles_acc_pPos_mc")->Clone();
+ TH1D* nPlutoElSts = (TH1D*) H1("fh_nof_particles_acc_pEl_sts")->Clone();
+ TH1D* nPlutoPosSts = (TH1D*) H1("fh_nof_particles_acc_pPos_sts")->Clone();
+ TH1D* nPlutoElRich = (TH1D*) H1("fh_nof_particles_acc_pEl_rich")->Clone();
+ TH1D* nPlutoPosRich = (TH1D*) H1("fh_nof_particles_acc_pPos_rich")->Clone();
+ TH1D* nPlutoElTrd = (TH1D*) H1("fh_nof_particles_acc_pEl_trd")->Clone();
+ TH1D* nPlutoPosTrd = (TH1D*) H1("fh_nof_particles_acc_pPos_trd")->Clone();
+ TH1D* nPlutoElTof = (TH1D*) H1("fh_nof_particles_acc_pEl_tof")->Clone();
+ TH1D* nPlutoPosTof = (TH1D*) H1("fh_nof_particles_acc_pPos_tof")->Clone();
+ TH1D* nUrqmdElMc2 = (TH1D*) H1("fh_nof_particles_acc_uEl_mc")->Clone();
+ TH1D* nUrqmdPosMc2 = (TH1D*) H1("fh_nof_particles_acc_uPos_mc")->Clone();
+ TH1D* nUrqmdElSts = (TH1D*) H1("fh_nof_particles_acc_uEl_sts")->Clone();
+ TH1D* nUrqmdPosSts = (TH1D*) H1("fh_nof_particles_acc_uPos_sts")->Clone();
+ TH1D* nUrqmdElRich = (TH1D*) H1("fh_nof_particles_acc_uEl_rich")->Clone();
+ TH1D* nUrqmdPosRich = (TH1D*) H1("fh_nof_particles_acc_uPos_rich")->Clone();
+ TH1D* nUrqmdElTrd = (TH1D*) H1("fh_nof_particles_acc_uEl_trd")->Clone();
+ TH1D* nUrqmdPosTrd = (TH1D*) H1("fh_nof_particles_acc_uPos_trd")->Clone();
+ TH1D* nUrqmdElTof = (TH1D*) H1("fh_nof_particles_acc_uEl_tof")->Clone();
+ TH1D* nUrqmdPosTof = (TH1D*) H1("fh_nof_particles_acc_uPos_tof")->Clone();
+
+ int n = 1;
+ int nRebinAcc = 4;
+ nPlutoElMc2->Rebin(n * nRebinAcc);
+ nPlutoElMc2->Scale(1. / (n * nRebinAcc));
+ nPlutoPosMc2->Rebin(n * nRebinAcc);
+ nPlutoPosMc2->Scale(1. / (n * nRebinAcc));
+ nPlutoElSts->Rebin(n * nRebinAcc);
+ nPlutoElSts->Scale(1. / (n * nRebinAcc));
+ nPlutoPosSts->Rebin(n * nRebinAcc);
+ nPlutoPosSts->Scale(1. / (n * nRebinAcc));
+ nPlutoElRich->Rebin(n * nRebinAcc);
+ nPlutoElRich->Scale(1. / (n * nRebinAcc));
+ nPlutoPosRich->Rebin(n * nRebinAcc);
+ nPlutoPosRich->Scale(1. / (n * nRebinAcc));
+ nPlutoElTrd->Rebin(n * nRebinAcc);
+ nPlutoElTrd->Scale(1. / (n * nRebinAcc));
+ nPlutoPosTrd->Rebin(n * nRebinAcc);
+ nPlutoPosTrd->Scale(1. / (n * nRebinAcc));
+ nPlutoElTof->Rebin(n * nRebinAcc);
+ nPlutoElTof->Scale(1. / (n * nRebinAcc));
+ nPlutoPosTof->Rebin(n * nRebinAcc);
+ nPlutoPosTof->Scale(1. / (n * nRebinAcc));
+
+ nUrqmdElMc2->Rebin(n * nRebinAcc);
+ nUrqmdElMc2->Scale(1. / (n * nRebinAcc));
+ nUrqmdPosMc2->Rebin(n * nRebinAcc);
+ nUrqmdPosMc2->Scale(1. / (n * nRebinAcc));
+ nUrqmdElSts->Rebin(n * nRebinAcc);
+ nUrqmdElSts->Scale(1. / (n * nRebinAcc));
+ nUrqmdPosSts->Rebin(n * nRebinAcc);
+ nUrqmdPosSts->Scale(1. / (n * nRebinAcc));
+ nUrqmdElRich->Rebin(n * nRebinAcc);
+ nUrqmdElRich->Scale(1. / (n * nRebinAcc));
+ nUrqmdPosRich->Rebin(n * nRebinAcc);
+ nUrqmdPosRich->Scale(1. / (n * nRebinAcc));
+ nUrqmdElTrd->Rebin(n * nRebinAcc);
+ nUrqmdElTrd->Scale(1. / (n * nRebinAcc));
+ nUrqmdPosTrd->Rebin(n * nRebinAcc);
+ nUrqmdPosTrd->Scale(1. / (n * nRebinAcc));
+ nUrqmdElTof->Rebin(n * nRebinAcc);
+ nUrqmdElTof->Scale(1. / (n * nRebinAcc));
+ nUrqmdPosTof->Rebin(n * nRebinAcc);
+ nUrqmdPosTof->Scale(1. / (n * nRebinAcc));
+
+ nPlutoElMc2->SetMinimum(min1);
+ nPlutoElMc2->SetMaximum(max1);
+ nPlutoPosMc2->SetMinimum(min1);
+ nPlutoPosMc2->SetMaximum(max1);
+ nPlutoElSts->SetMinimum(min1);
+ nPlutoElSts->SetMaximum(max1);
+ nPlutoPosSts->SetMinimum(min1);
+ nPlutoPosSts->SetMaximum(max1);
+ nPlutoElRich->SetMinimum(min1);
+ nPlutoElRich->SetMaximum(max1);
+ nPlutoPosRich->SetMinimum(min1);
+ nPlutoPosRich->SetMaximum(max1);
+ nPlutoElTrd->SetMinimum(min1);
+ nPlutoElTrd->SetMaximum(max1);
+ nPlutoPosTrd->SetMinimum(min1);
+ nPlutoPosTrd->SetMaximum(max1);
+ nPlutoElTof->SetMinimum(min1);
+ nPlutoElTof->SetMaximum(max1);
+ nPlutoPosTof->SetMinimum(min1);
+ nPlutoPosTof->SetMaximum(max1);
+
+ nUrqmdElMc2->SetMinimum(min1);
+ nUrqmdElMc2->SetMaximum(max1);
+ nUrqmdPosMc2->SetMinimum(min1);
+ nUrqmdPosMc2->SetMaximum(max1);
+ nUrqmdElSts->SetMinimum(min1);
+ nUrqmdElSts->SetMaximum(max1);
+ nUrqmdPosSts->SetMinimum(min1);
+ nUrqmdPosSts->SetMaximum(max1);
+ nUrqmdElRich->SetMinimum(min1);
+ nUrqmdElRich->SetMaximum(max1);
+ nUrqmdPosRich->SetMinimum(min1);
+ nUrqmdPosRich->SetMaximum(max1);
+ nUrqmdElTrd->SetMinimum(min1);
+ nUrqmdElTrd->SetMaximum(max1);
+ nUrqmdPosTrd->SetMinimum(min1);
+ nUrqmdPosTrd->SetMaximum(max1);
+ nUrqmdElTof->SetMinimum(min1);
+ nUrqmdElTof->SetMaximum(max1);
+ nUrqmdPosTof->SetMinimum(min1);
+ nUrqmdPosTof->SetMaximum(max1);
+
+ nPlutoElMc2->GetYaxis()->SetTitle("Yield");
+ nPlutoElMc2->SetTitle("Acceptance PLUTO particles");
+ nUrqmdElMc2->GetYaxis()->SetTitle("Yield");
+ nUrqmdElMc2->SetTitle("Acceptance UrQMD particles");
+
+ fHM->CreateCanvas("fh_nof_pluto_det", "fh_nof_pluto_det", 800, 800);
+ DrawH1({nPlutoElMc2, nPlutoPosMc2, nPlutoElSts, nPlutoPosSts, nPlutoElRich, nPlutoPosRich, nPlutoElTrd, nPlutoPosTrd,
+ nPlutoElTof, nPlutoPosTof},
+ {"", "", "", "", "", "", "", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofPP2 = new TLegend(0.45, 0.6, 0.88, 0.93);
+ legendNofPP2->SetFillColor(kWhite);
+ legendNofPP2->AddEntry(nPlutoElMc2, "electrons kMc");
+ legendNofPP2->AddEntry(nPlutoPosMc2, "positrons kMc");
+ legendNofPP2->AddEntry(nPlutoElSts, "electrons kAcc STS");
+ legendNofPP2->AddEntry(nPlutoPosSts, "positrons kAcc STS");
+ legendNofPP2->AddEntry(nPlutoElRich, "electrons kAcc STS + RICH");
+ legendNofPP2->AddEntry(nPlutoPosRich, "positrons kAcc STS + RICH");
+ legendNofPP2->AddEntry(nPlutoElTrd, "electrons kAcc STS + RICH + TRD");
+ legendNofPP2->AddEntry(nPlutoPosTrd, "positrons kAcc STS + RICH + TRD");
+ legendNofPP2->AddEntry(nPlutoElTof, "electrons kAcc STS + RICH + TRD + ToF");
+ legendNofPP2->AddEntry(nPlutoPosTof, "positrons kAcc STS + RICH + TRD + ToF");
+ legendNofPP2->Draw();
+
+ fHM->CreateCanvas("fh_nof_urqmd_det", "fh_nof_urqmd_det", 800, 800);
+ DrawH1({nUrqmdElMc2, nUrqmdPosMc2, nUrqmdElSts, nUrqmdPosSts, nUrqmdElRich, nUrqmdPosRich, nUrqmdElTrd, nUrqmdPosTrd,
+ nUrqmdElTof, nUrqmdPosTof},
+ {"", "", "", "", "", "", "", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofUP2 = new TLegend(0.45, 0.6, 0.88, 0.93);
+ legendNofUP2->SetFillColor(kWhite);
+ legendNofUP2->AddEntry(nUrqmdElMc2, "electrons kMc");
+ legendNofUP2->AddEntry(nUrqmdPosMc2, "positrons kMc");
+ legendNofUP2->AddEntry(nUrqmdElSts, "electrons kAcc STS");
+ legendNofUP2->AddEntry(nUrqmdPosSts, "positrons kAcc STS");
+ legendNofUP2->AddEntry(nUrqmdElRich, "electrons kAcc STS + RICH");
+ legendNofUP2->AddEntry(nUrqmdPosRich, "positrons kAcc STS + RICH");
+ legendNofUP2->AddEntry(nUrqmdElTrd, "electrons kAcc STS + RICH + TRD");
+ legendNofUP2->AddEntry(nUrqmdPosTrd, "positrons kAcc STS + RICH + TRD");
+ legendNofUP2->AddEntry(nUrqmdElTof, "electrons kAcc STS + RICH + TRD + ToF");
+ legendNofUP2->AddEntry(nUrqmdPosTof, "positrons kAcc STS + RICH + TRD + ToF");
+ legendNofUP2->Draw();
+
+ // Draw nofPoints of electrons and positrons in various detectors
+ TH1D* pElSts = (TH1D*) H1("fh_nof_points_pEl_sts")->Clone();
+ TH1D* pPosSts = (TH1D*) H1("fh_nof_points_pPos_sts")->Clone();
+ TH1D* pElRich = (TH1D*) H1("fh_nof_points_pEl_rich")->Clone();
+ TH1D* pPosRich = (TH1D*) H1("fh_nof_points_pPos_rich")->Clone();
+ TH1D* pElTrd = (TH1D*) H1("fh_nof_points_pEl_trd")->Clone();
+ TH1D* pPosTrd = (TH1D*) H1("fh_nof_points_pPos_trd")->Clone();
+ TH1D* pElTof = (TH1D*) H1("fh_nof_points_pEl_tof")->Clone();
+ TH1D* pPosTof = (TH1D*) H1("fh_nof_points_pPos_tof")->Clone();
+ TH1D* uElSts = (TH1D*) H1("fh_nof_points_uEl_sts")->Clone();
+ TH1D* uPosSts = (TH1D*) H1("fh_nof_points_uPos_sts")->Clone();
+ TH1D* uElRich = (TH1D*) H1("fh_nof_points_uEl_rich")->Clone();
+ TH1D* uPosRich = (TH1D*) H1("fh_nof_points_uPos_rich")->Clone();
+ TH1D* uElTrd = (TH1D*) H1("fh_nof_points_uEl_trd")->Clone();
+ TH1D* uPosTrd = (TH1D*) H1("fh_nof_points_uPos_trd")->Clone();
+ TH1D* uElTof = (TH1D*) H1("fh_nof_points_uEl_tof")->Clone();
+ TH1D* uPosTof = (TH1D*) H1("fh_nof_points_uPos_tof")->Clone();
+
+ pElSts->SetTitle("STS");
+ pElRich->SetTitle("RICH");
+ pElTrd->SetTitle("TRD");
+ pElTrd->GetXaxis()->SetRangeUser(0., 15.);
+ pElTof->SetTitle("ToF");
+ pElTof->GetXaxis()->SetRangeUser(0., 20.);
+
+ fHM->CreateCanvas("fh_nof_points_sts", "fh_nof_points_sts", 800, 800);
+ DrawH1({pElSts, pPosSts, uElSts, uPosSts}, {"", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofPSts = new TLegend(0.45, 0.7, 0.88, 0.93);
+ legendNofPSts->SetFillColor(kWhite);
+ legendNofPSts->AddEntry(pElSts, "PLUTO electrons");
+ legendNofPSts->AddEntry(pPosSts, "PLUTO positrons");
+ legendNofPSts->AddEntry(uElSts, "UrQMD electrons");
+ legendNofPSts->AddEntry(uPosSts, "UrQMD positrons");
+ legendNofPSts->Draw();
+
+ fHM->CreateCanvas("fh_nof_points_rich", "fh_nof_points_rich", 800, 800);
+ DrawH1({pElRich, pPosRich, uElRich, uPosRich}, {"", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99,
+ "hist p");
+
+ TLegend* legendNofPRich = new TLegend(0.45, 0.7, 0.88, 0.93);
+ legendNofPRich->SetFillColor(kWhite);
+ legendNofPRich->AddEntry(pElRich, "PLUTO electrons");
+ legendNofPRich->AddEntry(pPosRich, "PLUTO positrons");
+ legendNofPRich->AddEntry(uElRich, "UrQMD electrons");
+ legendNofPRich->AddEntry(uPosRich, "UrQMD positrons");
+ legendNofPRich->Draw();
+
+ fHM->CreateCanvas("fh_nof_points_trd", "fh_nof_points_trd", 800, 800);
+ DrawH1({pElTrd, pPosTrd, uElTrd, uPosTrd}, {"", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofPTrd = new TLegend(0.45, 0.7, 0.88, 0.93);
+ legendNofPTrd->SetFillColor(kWhite);
+ legendNofPTrd->AddEntry(pElTrd, "PLUTO electrons");
+ legendNofPTrd->AddEntry(pPosTrd, "PLUTO positrons");
+ legendNofPTrd->AddEntry(uElTrd, "UrQMD electrons");
+ legendNofPTrd->AddEntry(uPosTrd, "UrQMD positrons");
+ legendNofPTrd->Draw();
+
+ fHM->CreateCanvas("fh_nof_points_tof", "fh_nof_points_tof", 800, 800);
+ DrawH1({pElTof, pPosTof, uElTof, uPosTof}, {"", "", "", ""}, kLinear, kLog, false, 0.9, 0.8, 0.99, 0.99, "hist p");
+
+ TLegend* legendNofPTof = new TLegend(0.45, 0.7, 0.88, 0.93);
+ legendNofPTof->SetFillColor(kWhite);
+ legendNofPTof->AddEntry(pElTof, "PLUTO electrons");
+ legendNofPTof->AddEntry(pPosTof, "PLUTO positrons");
+ legendNofPTof->AddEntry(uElTof, "UrQMD electrons");
+ legendNofPTof->AddEntry(uPosTof, "UrQMD positrons");
+ legendNofPTof->Draw();
+
+ //Draw 2D Histos: Yield vs. Momentum and Pt
+ TH2D* pElMc = H2("fh_nof_plutoElectrons_p_pt_mc");
+ TH2D* pPosMc = H2("fh_nof_plutoPositrons_p_pt_mc");
+ TH2D* pElAcc = H2("fh_nof_plutoElectrons_p_pt_acc");
+ TH2D* pPosAcc = H2("fh_nof_plutoPositrons_p_pt_acc");
+
+ TH2D* uElMc = H2("fh_nof_urqmdElectrons_p_pt_mc");
+ TH2D* uPosMc = H2("fh_nof_urqmdPositrons_p_pt_mc");
+ TH2D* uElAcc = H2("fh_nof_urqmdElectrons_p_pt_acc");
+ TH2D* uPosAcc = H2("fh_nof_urqmdPositrons_p_pt_acc");
+
+ TCanvas* cPMc = fHM->CreateCanvas("fh_nof_plutoParticles_p_pt_mc", "fh_nof_plutoParticles_p_pt_mc", 1000, 500);
+ cPMc->Divide(2, 1);
+ cPMc->cd(1);
+ DrawH2(pElMc);
+ pElMc->SetMinimum(1e-6);
+ pElMc->SetMaximum(1e-1);
+ gPad->SetLogz(true);
+ cPMc->cd(2);
+ DrawH2(pPosMc);
+ pPosMc->SetMinimum(1e-6);
+ pPosMc->SetMaximum(1e-1);
+ gPad->SetLogz(true);
+
+ TCanvas* cPAcc = fHM->CreateCanvas("fh_nof_plutoParticles_p_pt_acc", "fh_nof_plutoParticles_p_pt_acc", 1000, 500);
+ cPAcc->Divide(2, 1);
+ cPAcc->cd(1);
+ DrawH2(pElAcc);
+ pElAcc->SetMinimum(1e-6);
+ pElAcc->SetMaximum(5e-2);
+ gPad->SetLogz(true);
+ cPAcc->cd(2);
+ DrawH2(pPosAcc);
+ pPosAcc->SetMinimum(1e-6);
+ pPosAcc->SetMaximum(5e-2);
+ gPad->SetLogz(true);
+
+ TCanvas* cUMc = fHM->CreateCanvas("fh_nof_urqmdParticles_p_pt_mc", "fh_nof_urqmdParticles_p_pt_mc", 1000, 500);
+ cUMc->Divide(2, 1);
+ cUMc->cd(1);
+ DrawH2(uElMc);
+ uElMc->SetMinimum(1e-6);
+ uElMc->SetMaximum(250);
+ gPad->SetLogz(true);
+ cUMc->cd(2);
+ DrawH2(uPosMc);
+ uPosMc->SetMinimum(1e-6);
+ uPosMc->SetMaximum(250);
+ gPad->SetLogz(true);
+
+ TCanvas* cUAcc = fHM->CreateCanvas("fh_nof_urqmdParticles_p_pt_acc", "fh_nof_urqmdParticles_p_pt_acc", 1000, 500);
+ cUAcc->Divide(2, 1);
+ cUAcc->cd(1);
+ DrawH2(uElAcc);
+ uElAcc->SetMinimum(1e-6);
+ uElAcc->SetMaximum(1e-1);
+ gPad->SetLogz(true);
+ cUAcc->cd(2);
+ DrawH2(uPosAcc);
+ uPosAcc->SetMinimum(1e-6);
+ uPosAcc->SetMaximum(1e-1);
+ gPad->SetLogz(true);
}
void CbmAnaDielectronTaskDraw::RebinMinvHist()
{
- int nRebin = 20;
+ int nRebin = 20; // general rebin
+ int nRebCB = 10; // rebin for CB histos
+ int nRebSP = 4; // rebin for single particle histos
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_combPairsPM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
+ H1("fh_combPairsPP_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
+ H1("fh_combPairsMM_minv_sameEvent_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
+ H1("fh_combPairsPM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
+ H1("fh_combPairsPP_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
+ H1("fh_combPairsMM_minv_mixedEvents_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebCB);
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);
@@ -263,14 +546,14 @@ 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);
+ H1("fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebSP);
+ H1("fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebSP);
+ H1("fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebSP);
+ H1("fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[i])->Rebin(nRebSP);
+ H1("fh_nof_plutoElectrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebSP);
+ H1("fh_nof_plutoPositrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebSP);
+ H1("fh_nof_urqmdElectrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebSP);
+ H1("fh_nof_urqmdPositrons_" + CbmLmvmHist::fAnaSteps[i])->Scale(1. / nRebSP);
for (int iP = 0; iP < CbmLmvmHist::fNofBgPairSources; iP++) {
stringstream ss;
@@ -390,8 +673,8 @@ void CbmAnaDielectronTaskDraw::SOverBg(CbmLmvmAnalysisSteps step)
{
TH1D* s = H1("fh_signal_minv_" + CbmLmvmHist::fAnaSteps[step]);
TH1D* bg = H1("fh_bg_minv_" + CbmLmvmHist::fAnaSteps[step]);
- H2("fh_signal_pty_" + CbmLmvmHist::fAnaSteps[step]);
- H2("fh_signal_pty_" + CbmLmvmHist::fAnaSteps[0]);
+ //H2("fh_signal_pty_" + CbmLmvmHist::fAnaSteps[step]); // TODO: commented these lines, what do they do??
+ //H2("fh_signal_pty_" + CbmLmvmHist::fAnaSteps[0]);
if (s->GetEntries() < 1) return;
diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
index 5b720d2571..9fe8d94700 100644
--- a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
+++ b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTaskDrawAll.cxx
@@ -80,12 +80,6 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
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);
@@ -94,7 +88,6 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
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);
@@ -104,7 +97,7 @@ void CbmAnaDielectronTaskDrawAll::DrawHistosFromFile(const string& fileNameInmed
fh_mean_nof_urqmdElectrons.resize(CbmLmvmHist::fNofAnaSteps);
fh_mean_nof_urqmdPositrons.resize(CbmLmvmHist::fNofAnaSteps);
- FillMeanHist();
+ FillMeanHist(); // TODO: Add method RebinHist() after this to rebin all histograms there to have better overview??
FillSumSignalsHist();
CalcCutEffRange(0.0, 0.2);
CalcCutEffRange(0.2, 0.6);
@@ -388,25 +381,27 @@ void CbmAnaDielectronTaskDrawAll::DrawMinv(CbmLmvmAnalysisSteps step)
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 *
- ******************************************************************************/
+ /*******************************************************************************
+ * 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;
+ bool setMinMax = true;
+
+ // Rebin factor for single particle and CB histograms. Their nBin in Task.cxx differ from each other
/* 1) draw number of e+ and e- per event */
// vs. momentum
@@ -435,57 +430,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
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);
@@ -620,12 +564,29 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
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();
+ TH1D* h21PMElidRaw = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone();
+ TH1D* h21PPElidRaw = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[kElId]->Clone();
+ TH1D* h21MMElidRaw = (TH1D*) fh_mean_combPairsMM_sameEvent_minv[kElId]->Clone();
+ TH1D* h21pmElidRaw = (TH1D*) fh_mean_combPairsPM_mixedEvents_minv[kElId]->Clone();
+ TH1D* h21ppElidRaw = (TH1D*) fh_mean_combPairsPP_mixedEvents_minv[kElId]->Clone();
+ TH1D* h21mmElidRaw = (TH1D*) fh_mean_combPairsMM_mixedEvents_minv[kElId]->Clone();
+
+ double bW = h21PMElidRaw->GetBinWidth(1); // MIND: 'bW' is used throughout this draw method!
+ cout << "DrawMinvCombSignalAndBg(): bW = " << bW << endl;
+
+ 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 << "DrawMinvCombSignalAndBg: nofEvents = " << nofEvents << endl;
+ }
+
+ h21PMElidRaw->Scale(bW * nofEvents);
+ h21PPElidRaw->Scale(bW * nofEvents);
+ h21MMElidRaw->Scale(bW * nofEvents);
+ h21pmElidRaw->Scale(bW * nofEvents);
+ h21ppElidRaw->Scale(bW * nofEvents);
+ h21mmElidRaw->Scale(bW * nofEvents);
h21PMElidRaw->GetXaxis()->SetTitle(xTitle.c_str());
h21PMElidRaw->GetYaxis()->SetTitle("absolute number");
@@ -659,8 +620,6 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
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());
@@ -817,9 +776,10 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
cocBg1Pt->Add(fh_mean_bg_minv[kPtCut]);
cocBg1Elid->Rebin(nRebin);
- cocBg1Elid->Scale(1. / bW);
+ double bWCoc = cocBg1Elid->GetBinWidth(1);
+ cocBg1Elid->Scale(1. / bWCoc);
cocBg1Pt->Rebin(nRebin);
- cocBg1Pt->Scale(1. / bW);
+ cocBg1Pt->Scale(1. / bWCoc);
Bpm1Elid->GetXaxis()->SetRangeUser(0, 2.);
Bpm1Elid->GetXaxis()->SetTitle(xTitle.c_str());
@@ -1203,31 +1163,35 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
TH1D* h63SigNpmTt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kTtCut]->Clone();
TH1D* h63SigNpmPt = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kPtCut]->Clone();
+ TH1D* h63SigNpmPt2 = (TH1D*) fh_mean_combSignalNpm_assemb_minv[kPtCut]->Clone();
+
if (setMinMax) {
- h63SigNpmElid->SetMinimum(yMin / 10.); // setting min/max somehow makes the histos empty
+ h63SigNpmElid->SetMinimum(yMin); // setting min/max somehow makes the histos empty
h63SigNpmElid->SetMaximum(yMax);
- h63SigNpmGamma->SetMinimum(yMin / 10.);
+ h63SigNpmGamma->SetMinimum(yMin);
h63SigNpmGamma->SetMaximum(yMax);
- h63SigNpmMvd1->SetMinimum(yMin / 10.);
+ h63SigNpmMvd1->SetMinimum(yMin);
h63SigNpmMvd1->SetMaximum(yMax);
- h63SigNpmSt->SetMinimum(yMin / 10.);
+ h63SigNpmSt->SetMinimum(yMin);
h63SigNpmSt->SetMaximum(yMax);
- h63SigNpmRt->SetMinimum(yMin / 10.);
+ h63SigNpmRt->SetMinimum(yMin);
h63SigNpmRt->SetMaximum(yMax);
- h63SigNpmTt->SetMinimum(yMin / 10.);
+ h63SigNpmTt->SetMinimum(yMin);
h63SigNpmTt->SetMaximum(yMax);
- h63SigNpmPt->SetMinimum(yMin / 10.);
+ h63SigNpmPt->SetMinimum(yMin);
h63SigNpmPt->SetMaximum(yMax);
+ h63SigNpmPt2->SetMinimum(yMin);
+ h63SigNpmPt2->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");*/
+ h63SigNpmPt2->GetXaxis()->SetRangeUser(0, 2.);
+ h63SigNpmPt2->GetXaxis()->SetTitle(xTitle.c_str());
+ h63SigNpmPt2->GetYaxis()->SetTitle(yTitle.c_str());
+ h63SigNpmPt2->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},
@@ -1267,13 +1231,14 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
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");
+ 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({h63SigNpmPt2}, {""}, 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();*/
+ legend63e4->AddEntry(h63SigNpmPt2, "P_{t} cut");
+ legend63e4->Draw();
// from 'Cocktail + BG'
// elid
@@ -1303,24 +1268,27 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
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");
+ 0.99, 0.99, "");
+
+ /*h65coc1Elid->SetFillColor(kMagenta - 3);
+ h65coc1Elid->SetLineColor(kMagenta - 2);
+ h65coc1Elid->SetLineStyle(1);
+ h65coc1Elid->SetLineWidth(3);
+ h65coc1Elid->SetFillStyle(3344);*/
+
+ h65cocBgElid->SetLineStyle(0);
+ h65cbAssElid->SetLineStyle(2);
+ h65SigCocBgElid->SetLineStyle(3);
+ h65coc1Elid->SetLineStyle(1);
TLegend* legend65 = new TLegend(0.55, 0.75, 0.98, 0.9);
- legend65->SetFillColor(kWhite);
+ //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();
@@ -1350,8 +1318,7 @@ void CbmAnaDielectronTaskDrawAll::DrawMinvCombSignalAndBg()
legend66->Draw();
/* 7) Error */
-
- TH1D* h71NpmElid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv_raw[kElId]->Clone();
+ TH1D* h71NpmElid = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[kElId]->Clone();
TH1D* h71sigElid = (TH1D*) fh_mean_combSignal_errProp_minv[kElId]->Clone();
h71NpmElid->GetXaxis()->SetRangeUser(0, 2.);
@@ -1517,13 +1484,6 @@ void CbmAnaDielectronTaskDrawAll::FillMeanHist()
}
}
- 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;
@@ -1621,8 +1581,7 @@ void CbmAnaDielectronTaskDrawAll::CalcCutEffRange(Double_t minMinv, Double_t max
void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
{
-
- int nofBinsRaw;
+ int nofBins;
double bW = 0;
Int_t nofEvents = 0;
for (int i = 0; i < fNofSignals; i++) {
@@ -1648,14 +1607,16 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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_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);
+ fh_mean_combPairsMM_mixedEvents_minv_raw[step]->Rebin(nRebin);*/
- nofBinsRaw = fh_mean_combPairsPM_sameEvent_minv_raw[step]->GetNbinsX();
+ //int nofBins = hBpp->GetNbinsX(); TODO: test; delete
+ //nofBinsRaw = fh_mean_combPairsPM_sameEvent_minv_raw[step]->GetNbinsX(); _RAW_
+ nofBins = fh_mean_combPairsPM_sameEvent_minv[step]->GetNbinsX(); // MIND: 'nofBins' is used throughout this method!
// calculate geom. mean of same events
TH1D* hBpp = (TH1D*) fh_mean_combPairsPP_sameEvent_minv[step]->Clone();
@@ -1663,8 +1624,6 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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);
@@ -1714,9 +1673,9 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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;
+ double k2 = fh_mean_combBg_k_minv[step]->GetBinContent(iBin);
+ double gm2 = fh_mean_combBg_GeomMeanMixed_minv[step]->GetBinContent(iBin);
+ double content = 2 * k2 * gm2 * normGM;
fh_mean_combBg_assemb_minv[step]->SetBinContent(iBin, content);
}
@@ -1725,9 +1684,9 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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'
+ TH1D* cock;
if (step == kMc) cock = (TH1D*) GetCoctailMinv(kMc);
else if (step == kAcc)
cock = (TH1D*) GetCoctailMinv(kAcc);
@@ -1771,18 +1730,18 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
fh_mean_combSignalBCoc_assemb_minv[step]->Add(cbAss, -1.);
// error calculation
- TH1D* err = (TH1D*) fh_mean_combPairsPM_sameEvent_minv_raw[step]->Clone();
+ TH1D* err = (TH1D*) fh_mean_combPairsPM_sameEvent_minv[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);
+ for (int iBin = 1; iBin <= nofBins; iBin++) { // calculate error propagation
+ double Npm = fh_mean_combPairsPM_sameEvent_minv[step]->GetBinContent(iBin);
+ double Bpp = fh_mean_combPairsPP_sameEvent_minv[step]->GetBinContent(iBin);
+ double Bmm = fh_mean_combPairsMM_sameEvent_minv[step]->GetBinContent(iBin);
+ double bpm = fh_mean_combPairsPM_mixedEvents_minv[step]->GetBinContent(iBin);
+ double bpp = fh_mean_combPairsPP_mixedEvents_minv[step]->GetBinContent(iBin);
+ double bmm = fh_mean_combPairsMM_mixedEvents_minv[step]->GetBinContent(iBin);
// calculation of error propagation
double DNpm = TMath::Sqrt(Npm); // Δ<B+->
@@ -1829,15 +1788,15 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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)) {
+ if (iBin == 3 || iBin == 4 || ((iBin % (500 / nRebin) == 0 && iBin <= (2000 / nRebin)))) {
cout << "step = " << step << endl;
+ cout << "iBin = " << iBin << 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;
@@ -1859,12 +1818,13 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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;
+ cout << "Bin error Npm = " << fh_mean_combPairsPM_sameEvent_minv[step]->GetBinError(iBin) << endl;
}
} // error propagation
// scale histograms
bW = fh_mean_combPairsPM_sameEvent_minv[step]->GetBinWidth(1);
+ cout << "bW = " << bW << endl;
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));
@@ -1877,9 +1837,23 @@ void CbmAnaDielectronTaskDrawAll::CalcCombBGHistos()
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
+
+ for (int iBin = 1; iBin <= nofBins; iBin++) { // TODO: only to check error value after sacling
+ if (iBin == 3 || iBin == 4 || (iBin % (500 / nRebin) == 0 && iBin <= (2000 / nRebin))) {
+ double Npm = fh_mean_combPairsPM_sameEvent_minv[step]->GetBinContent(iBin);
+ double Bm = fh_mean_combPairsMM_sameEvent_minv[step]->GetBinContent(iBin);
+ cout << "CHECK: step = " << step << endl;
+ cout << "CHECK: iBin = " << iBin << endl;
+ cout << "CHECK: Npm = " << Npm * nofEvents * bW << endl;
+ cout << "CHECK: Bm = " << Bm * nofEvents * bW << endl;
+ cout << "CHECK: Bin error Npm = " << fh_mean_combPairsPM_sameEvent_minv[step]->GetBinError(iBin)
+ << endl;
+ cout << "CHECK: Bin error Sig_ass (Npm) = " << fh_mean_combSignalNpm_assemb_minv[step]->GetBinError(iBin)
+ << endl;
+ }
+ }
+ } // steps
}
/*void CbmAnaDielectronTaskDrawAll::CompareSTSversions()
@@ -1991,7 +1965,7 @@ void CbmAnaDielectronTaskDrawAll::DrawSBgSignals()
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* bg = (TH1D*) fh_mean_bg_minv[step]->Clone(); // TODO: not better to take BG of PP intead of mean BG??
Int_t nofEvents = (int) H1(iF, "fh_event_number")->GetEntries();
s->Scale(1. / nofEvents);
cFit->cd();
diff --git a/macro/analysis/dielectron/batch_job.py b/macro/analysis/dielectron/batch_job.py
old mode 100644
new mode 100755
index f76179266b..b61b1c4caa
--- a/macro/analysis/dielectron/batch_job.py
+++ b/macro/analysis/dielectron/batch_job.py
@@ -1,38 +1,70 @@
-#!/ 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 = 100
- 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')
+ jobId = os.environ.get('SLURM_ARRAY_JOB_ID')
- taskId = os.environ.get('SLURM_ARRAY_TASK_ID')
+ colEnergy = "8gev"
+ colSystem = "auau"
- jobId = os.environ.get('SLURM_ARRAY_JOB_ID') colEnergy = "8gev" colSystem = "auau" print("dataDir:" + dataDir)
+ print("dataDir:" + dataDir)
+ os.system(("source {}").format(cbmrootConfigPath))
- os.system(("source {}").format(cbmrootConfigPath))
+ workDir = dataDir + "/workdir/" + jobId + "_" + taskId + "/"
+ if os.path.exists(workDir):
+ shutil.rmtree(workDir)
+ os.makedirs(workDir)
+ os.chdir(workDir)
- 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)
-#plutoFile = ""
- plutoFile = getPlutoPath(colSystem, colEnergy, plutoParticle, taskId) #this one was activated before
-#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"
- 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"
+ print("Here: in batchJobPy")
+ #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))
- 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))
+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 ""
-
- 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()
+if __name__ == '__main__':
+ main()
diff --git a/macro/analysis/dielectron/batch_job_common.py b/macro/analysis/dielectron/batch_job_common.py
old mode 100644
new mode 100755
index 2e49d9406f..90b06b28e2
--- a/macro/analysis/dielectron/batch_job_common.py
+++ b/macro/analysis/dielectron/batch_job_common.py
@@ -1,58 +1,94 @@
-#!/ 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 ""
-
+#!/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"
+ colEnergyDir = "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 == "rho0":
+ offset = 4000
+ 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 do not contain '/TGeant4/'
+ if plutoParticle == "w" or plutoParticle == "wdalitz":
+ mcFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".tra.root"
+ parFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".par.root"
+ digiFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".event.raw.root"
+ recoFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/" + str(offsetId) + "/" + str(offsetId) + ".rec.root"
+ else:
+ mcFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".tra.root"
+ parFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".par.root"
+ digiFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".event.raw.root"
+ recoFile = recoDataDirIn + colSystem + "/" + colEnergyDir + "/centr_0_10/sis100_electron_target_25_mkm/TGeant4/" + str(offsetId) + "/" + str(offsetId) + ".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 100644
new mode 100755
index 60223be59b..76e9230106
--- a/macro/analysis/dielectron/batch_send.py
+++ b/macro/analysis/dielectron/batch_send.py
@@ -1,29 +1,53 @@
-#!/ 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()
+#!/usr/bin/env python3
+
+import os
+import shutil
+
+def main():
+ nofJobs = 100
+ 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
old mode 100644
new mode 100755
index 288932097a..a76bca4e44
--- a/macro/analysis/dielectron/batch_send_common.py
+++ b/macro/analysis/dielectron/batch_send_common.py
@@ -1,4 +1,4 @@
-#!/ usr / bin / env python3
+#!/usr/bin/env python3
## #changed on Feb 8th; first check if works or not ###
@@ -9,14 +9,14 @@ def main():
nofJobs = 100
timeLimit="08:00:00"
geoSetup="sis100_electron"
- plutoParticles = ["inmed_had_epem", "wdalitz", "w", "phi", "qgp_epem"]
+ plutoParticles = ["rho0", "wdalitz", "w", "phi", "qgp_epem"]
- recoDataDirIn = "/lustre/cbm/pwg/common/mc/cbmsim/apr20_fr_18.2.1_fs_jun19p1/urqmd_pluto_"
+ recoDataDirIn = "/lustre/nyx/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
+ #All data in data dir will be removed
removeData = False
if removeData:
print("All data in dataDirOut will be removed. Dir:" + dataDirOut)
@@ -27,13 +27,13 @@ def main():
for plutoParticle in plutoParticles:
recoDataDirInPluto = recoDataDirIn + plutoParticle + "/"
- dataDirOutPluto = dataDirOut + 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
+ #- 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)
diff --git a/macro/analysis/dielectron/draw_all.py b/macro/analysis/dielectron/draw_all.py
old mode 100644
new mode 100755
index 2de2035e4b..aa6a559486
--- a/macro/analysis/dielectron/draw_all.py
+++ b/macro/analysis/dielectron/draw_all.py
@@ -1,20 +1,29 @@
-#!/ usr / bin / env python3
+#!/usr/bin/env python3
import os, shutil
- 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()
+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_analysis.C b/macro/analysis/dielectron/draw_analysis.C
old mode 100644
new mode 100755
diff --git a/macro/analysis/dielectron/draw_analysis_all.C b/macro/analysis/dielectron/draw_analysis_all.C
old mode 100644
new mode 100755
diff --git a/macro/analysis/dielectron/hadd_many.py b/macro/analysis/dielectron/hadd_many.py
old mode 100644
new mode 100755
index 249f3357d1..1ac5f91805
--- a/macro/analysis/dielectron/hadd_many.py
+++ b/macro/analysis/dielectron/hadd_many.py
@@ -1,11 +1,24 @@
-import os import shutil
+import os
+import shutil
- def main() : plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"] dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm"
-
- for plutoParticle in plutoParticles : dataDirPluto = dataDir + "/" + plutoParticle
-
- 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))
-
- 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()
+def main():
+ plutoParticles =["inmed", "omegadalitz", "omegaepem", "phi", "qgp"]
+ dataDir = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm"
+
+ for plutoParticle in plutoParticles:
+ dataDirPluto = dataDir + "/" + plutoParticle
+
+ #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))
+
+ 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_analysis.C b/macro/analysis/dielectron/run_analysis.C
old mode 100644
new mode 100755
index a387f76bbf..d5718321c7
--- a/macro/analysis/dielectron/run_analysis.C
+++ b/macro/analysis/dielectron/run_analysis.C
@@ -8,8 +8,20 @@ void run_analysis(const string& mcFile = "/lustre/nyx/cbm/users/criesen/c
const string& recoFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/reco.1.root",
const string& analysisFile = "/lustre/nyx/cbm/users/criesen/cbm/data/lmvm/inmed/analysis.1.root",
const string& plutoParticle = "inmed", const string& colSystem = "auau",
- const string& colEnergy = "8gev", const string& geoSetup = "sis100_electron", int nEvents = 100)
+ const string& colEnergy = "8gev", const string& geoSetup = "sis100_electron", int nEvents = 1000)
{
+ /*cout << "Here: just entered run_ana.C" << endl;
+ Int_t Interval=10;
+ Int_t PID=gSystem->GetPid();
+ cout<<"PID: "<<PID<<endl;
+ TString cmdline="$VMCWORKDIR/macro/analysis/dielectron/check_memory.sh ";
+ cmdline+= PID;
+ cmdline+= " ";
+ cmdline+= Interval;
+ cmdline+= " &";
+ cout<<cmdline<<endl;
+ gSystem->Exec(cmdline);*/
+
TTree::SetMaxTreeSize(90000000000);
@@ -61,7 +73,7 @@ void run_analysis(const string& mcFile = "/lustre/nyx/cbm/users/criesen/c
FairRuntimeDb* rtdb = run->GetRuntimeDb();
FairParRootFileIo* parIo1 = new FairParRootFileIo();
FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
- parIo1->open(parFile.c_str(), "UPDATE");
+ parIo1->open(parFile.c_str(), "READ");
rtdb->setFirstInput(parIo1);
if (!parFileList->IsEmpty()) {
parIo2->open(parFileList, "in");
@@ -69,14 +81,12 @@ void run_analysis(const string& mcFile = "/lustre/nyx/cbm/users/criesen/c
}
std::cout << std::endl << "-I- " << myName << ": Initialise run" << std::endl;
run->Init();
-
rtdb->setOutput(parIo1);
rtdb->saveOutput();
rtdb->print();
std::cout << "-I- " << myName << ": Starting run" << std::endl;
run->Run(0, nEvents);
-
timer.Stop();
std::cout << std::endl << std::endl;
std::cout << "Macro finished succesfully." << std::endl;
diff --git a/macro/analysis/dielectron/run_litqa.C b/macro/analysis/dielectron/run_litqa.C
index 8808346e0f..7a3020d00b 100644
--- a/macro/analysis/dielectron/run_litqa.C
+++ b/macro/analysis/dielectron/run_litqa.C
@@ -112,7 +112,7 @@ void run_litqa(const string& mcFile = "/lustre/nyx/cbm/users/criesen/cbm/data/
FairRuntimeDb* rtdb = run->GetRuntimeDb();
FairParRootFileIo* parIo1 = new FairParRootFileIo();
FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
- parIo1->open(parFile.c_str(), "UPDATE");
+ parIo1->open(parFile.c_str(), "READ");
rtdb->setFirstInput(parIo1);
if (!parFileList->IsEmpty()) {
parIo2->open(parFileList, "in");
--
GitLab