diff --git a/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx b/analysis/PWGDIL/dielectron/lmvm/CbmAnaDielectronTask.cxx index e5fb3add05a44a4ff47345b5ade8ec7e8dd452be..677a3686e2e818b201a81d162873732ecd6fc9e7 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 d46ec889f487e405cbc85a29770979ddbc788635..7926e15075c15db6732873e50143f2ddd70ba4be 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 c3d3b24ad4a82f406a6e49a28b2b14ded0681659..119013da6112e11eb8fface5f807a1e18d225a34 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 5b720d2571c439a8ee34059a16812755e6d29b0d..9fe8d947008e6e18ae73efc7e15fbf1305c6fea8 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 f76179266b26a1f0b9f15b45778e351459acecfd..b61b1c4caa360cfac340969a66a92e11f409af5b --- 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 2e49d9406f9dd0a57235de2f75e501cc5958a708..90b06b28e2c350ffad16d017035d640206b6aaae --- 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 60223be59bae3a4eb99a198002522aefd22754c6..76e92301068d406eb209a6267f5bb6583f7be93b --- 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 288932097ac2315195ed7f8cb4c267c178069c69..a76bca4e44763dd3e989cabee150f94b3287151b --- 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 2de2035e4b349df2b192af31f6be8c431bf9dad2..aa6a559486f46848f7067754b865d0fef0919078 --- 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 249f3357d1f4045141d9e76559942aca933bc70f..1ac5f9180558a3742950b12fee93b802e5bbbd9c --- 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 a387f76bbfd6211310df4f10057c8b652a9ac73f..d5718321c713802bf06083f3a4516e9a2b89a690 --- 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 8808346e0f351c2685da3515ac3a0d597ffb1364..7a3020d00bca7caa77bc03778ddebec6d2313b6e 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");