diff --git a/PWGEM/Dilepton/Tasks/lmeeLFCocktail.cxx b/PWGEM/Dilepton/Tasks/lmeeLFCocktail.cxx index ca760249fd4..57f6761aa65 100644 --- a/PWGEM/Dilepton/Tasks/lmeeLFCocktail.cxx +++ b/PWGEM/Dilepton/Tasks/lmeeLFCocktail.cxx @@ -100,11 +100,6 @@ struct lmeelfcocktail { TH1F* fhwMultmT2; TH1F* fhKW; TF1* ffVPHpT; - TObjArray* fArr; - TObjArray* fArrResoPt; - TObjArray* fArrResoEta; - TObjArray* fArrResoPhi_Pos; - TObjArray* fArrResoPhi_Neg; std::vector> fmee_orig, fmotherpT_orig, fphi_orig, frap_orig, fmee_orig_wALT, fmotherpT_orig_wALT, fmee, fphi, frap, fmee_wALT; std::vector> fpteevsmee_wALT, fpteevsmee_orig_wALT, fpteevsmee_orig, fpteevsmee; @@ -135,14 +130,9 @@ struct lmeelfcocktail { Configurable fConfigResFileName{"cfgResFileName", "", "name of resolution file"}; Configurable fConfigEffFileName{"cfgEffFileName", "", "name of efficiency file"}; Configurable fConfigMinOpAng{"cfgMinOpAng", 0.050, "minimum opening angle"}; - Configurable fConfigDoRapidityCut{"cfgDoRapidityCut", false, "apply rapidity cut"}; - Configurable fConfigRapidityCut{"cfgRapidityCut", 1.0, "rapdity cut"}; Configurable fConfigNBinsPhi{"cfgNBinsPhi", 240, "number of bins in phi"}; - Configurable fConfigMinPhi{"cfgMinPhi", 0.0, "lowerst bin in phi"}; - Configurable fConfigMaxPhi{"cfgMaxPhi", TMath::TwoPi(), "hightest bin in phi"}; Configurable fConfigNBinsRap{"cfgNBinsRap", 240, "number of bins in rap"}; - Configurable fConfigMinRap{"cfgMinRap", -1.2, "lowest bin in rap"}; - Configurable fConfigMaxRap{"cfgMaxRap", 1.2, "hightest bin in rap"}; + Configurable fConfigMaxAbsRap{"cfgMaxAbsRap", 1.2, "bin range in rap"}; Configurable fConfigEffHistName{"cfgEffHistName", "fhwEffpT", "hisogram name in efficiency file"}; Configurable fConfigResPHistName{"cfgResPHistName", "ptSlices", "histogram name for p in resolution file"}; Configurable fConfigResPtHistName{"cfgResPtHistName", "RelPtResArrCocktail", "histogram name for pt in resolution file"}; @@ -156,17 +146,15 @@ struct lmeelfcocktail { Configurable fConfigMultHistPt2Name{"cfgMultHistPt2Name", "fhwMultpT_upperlimit", "histogram name for pt 2 in multiplicity file"}; Configurable fConfigMultHistMtName{"cfgMultHistMtName", "fhwMultmT", "histogram name for mt in multiplicity file"}; Configurable fConfigMultHistMt2Name{"cfgMultHistMt2Name", "fhwMultmT_upperlimit", "histogram name for mt 2 in multiplicity file"}; - Configurable fConfigKWNBins{"cfgKWNBins", 10000, "number of bins for Kroll-Wada"}; Configurable fConfigKWMax{"cfgKWMax", 1.1, "upper bound of Kroll-Wada"}; Configurable fConfigDoVirtPh{"cfgDoVirtPh", false, "generate one virt. photon for each pion"}; Configurable fConfigPhotonPtFileName{"cfgPhotonPtFileName", "", "file name for photon pT parametrization"}; Configurable fConfigPhotonPtDirName{"cfgPhotonPtDirName", "", "directory name for photon pT parametrization"}; Configurable fConfigPhotonPtFuncName{"cfgPhotonPtFuncName", "111_pt", "function name for photon pT parametrization"}; - // Configurable axes crashed the task. Take them out for the moment - // ConfigurableAxis fConfigPtBins{"cfgPtBins", {VARIABLE_WIDTH, 0., 0.5, 1, 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5, 7., 7.5, 8.}, "pT bins"}; - // ConfigurableAxis fConfigMBins{"cfgMBins", {VARIABLE_WIDTH, 0., 0.08, 0.14, 0.2, 1.1, 2.7, 2.8, 3.2, 5.0}, "mee bins"}; - // ConfigurableAxis fConfigDCABins{"cfgDCABins", {VARIABLE_WIDTH, 0., 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 3., 4., 5., 7., 10.}, "DCA bins"}; + ConfigurableAxis fConfigPtBins{"cfgPtBins", {VARIABLE_WIDTH, 0., 0.5, 1, 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5, 7., 7.5, 8.}, "pT bins"}; + ConfigurableAxis fConfigMBins{"cfgMBins", {VARIABLE_WIDTH, 0., 0.08, 0.14, 0.2, 1.1, 2.7, 2.8, 3.2, 5.0}, "mee bins"}; + ConfigurableAxis fConfigDCABins{"cfgDCABins", {VARIABLE_WIDTH, 0., 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 3., 4., 5., 7., 10.}, "DCA bins"}; Configurable> fConfigDCATemplateEdges{"cfgDCATemplateEdges", {0., .3, .4, .6, 1., 2.}, "DCA template edges"}; @@ -186,9 +174,7 @@ struct lmeelfcocktail { } GetEffHisto(TString(fConfigEffFileName), TString(fConfigEffHistName)); - InitSmearer(TString(fConfigResFileName), TString(fConfigResPtHistName), TString(fConfigResEtaHistName), TString(fConfigResPhiPosHistName), TString(fConfigResPhiNegHistName)); - GetDCATemplates(TString(fConfigDCAFileName), TString(fConfigDCAHistName)); GetMultHisto(TString(fConfigMultFileName), TString(fConfigMultHistPtName), TString(fConfigMultHistPt2Name), TString(fConfigMultHistMtName), TString(fConfigMultHistMt2Name)); if (fConfigDoVirtPh) { @@ -202,500 +188,506 @@ struct lmeelfcocktail { void run(o2::framework::ProcessingContext& pc) { - // get the tracks - auto mctracks = pc.inputs().get>("mctracks"); - registry.fill(HIST("NEvents"), 0.5); - - std::vector eBuff; - std::vector echBuff; - std::vector eweightBuff; - - bool skipNext = false; - - int trackID = -1; - // Loop over all MC particle - for (auto& mctrack : mctracks) { - trackID++; - if (o2::mcgenstatus::getHepMCStatusCode(mctrack.getStatusCode()) != 1) - continue; - if (abs(mctrack.GetPdgCode()) == 11) { - // get the electron - //--------------- - if (fConfigDoPairing) { - // LS and ULS spectra - PxPyPzEVector e, dielectron; - Char_t ech, dielectron_ch; - Double_t eweight, dielectron_weight; - e.SetPxPyPzE(mctrack.Px(), mctrack.Py(), mctrack.Pz(), - mctrack.GetEnergy()); - if (mctrack.GetPdgCode() > 0) { - ech = 1.; - } else { - ech = -1.; - } - eweight = mctrack.getWeight(); - // put in the buffer - //----------------- - eBuff.push_back(e); - echBuff.push_back(ech); - eweightBuff.push_back(eweight); - // loop the buffer and pair - //------------------------ - for (Int_t jj = eBuff.size() - 2; jj >= 0; jj--) { - dielectron = eBuff.at(jj) + e; - dielectron_ch = (echBuff.at(jj) + ech) / 2; - dielectron_weight = eweightBuff.at(jj) * eweight; - - if (dielectron_ch == 0) - registry.fill(HIST("ULS_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); - if (dielectron_ch > 0) - registry.fill(HIST("LSpp_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); - if (dielectron_ch < 0) - registry.fill(HIST("LSmm_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); - if (e.Pt() > fConfigMinPt && eBuff.at(jj).Pt() > fConfigMinPt && e.Pt() < fConfigMaxPt && eBuff.at(jj).Pt() < fConfigMaxPt && TMath::Abs(e.Eta()) < fConfigMaxEta && TMath::Abs(eBuff.at(jj).Eta()) < fConfigMaxEta && e.Vect().Unit().Dot(eBuff.at(jj).Vect().Unit()) < TMath::Cos(fConfigMinOpAng)) { + // get number of events per timeframe + auto Nparts = pc.inputs().getNofParts(0); + + for (auto i = 0U; i < Nparts; ++i) { + registry.fill(HIST("NEvents"), 0.5); + // get the tracks + auto mctracks = pc.inputs().get>("mctracks", i); + + std::vector eBuff; + std::vector echBuff; + std::vector eweightBuff; + + bool skipNext = false; + + int trackID = -1; + // Loop over all MC particle + for (auto& mctrack : mctracks) { + trackID++; + if (o2::mcgenstatus::getHepMCStatusCode(mctrack.getStatusCode()) != 1) + continue; + if (abs(mctrack.GetPdgCode()) == 11) { + // get the electron + //--------------- + if (fConfigDoPairing) { + // LS and ULS spectra + PxPyPzEVector e, dielectron; + Char_t ech, dielectron_ch; + Double_t eweight, dielectron_weight; + e.SetPxPyPzE(mctrack.Px(), mctrack.Py(), mctrack.Pz(), + mctrack.GetEnergy()); + if (mctrack.GetPdgCode() > 0) { + ech = 1.; + } else { + ech = -1.; + } + eweight = mctrack.getWeight(); + // put in the buffer + //----------------- + eBuff.push_back(e); + echBuff.push_back(ech); + eweightBuff.push_back(eweight); + // loop the buffer and pair + //------------------------ + for (Int_t jj = eBuff.size() - 2; jj >= 0; jj--) { + dielectron = eBuff.at(jj) + e; + dielectron_ch = (echBuff.at(jj) + ech) / 2; + dielectron_weight = eweightBuff.at(jj) * eweight; + if (dielectron_ch == 0) - registry.fill(HIST("ULS"), dielectron.M(), dielectron.Pt(), dielectron_weight); + registry.fill(HIST("ULS_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); if (dielectron_ch > 0) - registry.fill(HIST("LSpp"), dielectron.M(), dielectron.Pt(), dielectron_weight); + registry.fill(HIST("LSpp_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); if (dielectron_ch < 0) - registry.fill(HIST("LSmm"), dielectron.M(), dielectron.Pt(), dielectron_weight); + registry.fill(HIST("LSmm_orig"), dielectron.M(), dielectron.Pt(), dielectron_weight); + if (e.Pt() > fConfigMinPt && eBuff.at(jj).Pt() > fConfigMinPt && e.Pt() < fConfigMaxPt && eBuff.at(jj).Pt() < fConfigMaxPt && TMath::Abs(e.Eta()) < fConfigMaxEta && TMath::Abs(eBuff.at(jj).Eta()) < fConfigMaxEta && e.Vect().Unit().Dot(eBuff.at(jj).Vect().Unit()) < TMath::Cos(fConfigMinOpAng)) { + if (dielectron_ch == 0) + registry.fill(HIST("ULS"), dielectron.M(), dielectron.Pt(), dielectron_weight); + if (dielectron_ch > 0) + registry.fill(HIST("LSpp"), dielectron.M(), dielectron.Pt(), dielectron_weight); + if (dielectron_ch < 0) + registry.fill(HIST("LSmm"), dielectron.M(), dielectron.Pt(), dielectron_weight); + } } } - } - - if (skipNext) { - skipNext = false; - continue; // skip if marked as second electron - } - if (!(mctrack.getMotherTrackId() > -1)) - continue; // has no mother - - auto const& mother = mctracks[mctrack.getMotherTrackId()]; - - if (mother.getMotherTrackId() > -1) - continue; // mother is not primary - - if (mctrack.getSecondMotherTrackId() - mctrack.getMotherTrackId() > 0) - continue; // more than one mother - - // skip for the moment other particles rather than pi0, eta, etaprime, - // omega, rho, phi. - switch (mother.GetPdgCode()) { - case 111: - break; - case 221: - break; - case 331: - break; - case 113: - break; - case 223: - break; - case 333: - break; - case 443: - break; - default: - continue; - } + if (skipNext) { + skipNext = false; + continue; // skip if marked as second electron + } - // Not sure about this cut. From GammaConv group. Harmless a priori. - if (!(fabs(mctrack.GetEnergy() - mctrack.Pz()) > 0.)) - continue; + if (!(mctrack.getMotherTrackId() > -1)) + continue; // has no mother + + auto const& mother = mctracks[mctrack.getMotherTrackId()]; + + if (mother.getMotherTrackId() > -1) + continue; // mother is not primary + + if (mctrack.getSecondMotherTrackId() - mctrack.getMotherTrackId() > 0) + continue; // more than one mother + + // skip for the moment other particles rather than pi0, eta, etaprime, + // omega, rho, phi. + switch (mother.GetPdgCode()) { + case 111: + break; + case 221: + break; + case 331: + break; + case 113: + break; + case 223: + break; + case 333: + break; + case 443: + break; + default: + continue; + } - // ???? this applied only to first daughter! - Double_t yPre = (mctrack.GetEnergy() + mctrack.Pz()) / (mctrack.GetEnergy() - mctrack.Pz()); - Double_t y = 0.5 * TMath::Log(yPre); - if (fConfigDoRapidityCut) { // Apply rapidity cut on mother consistent with GammaConv group. (??? but it is not applied on mother?) - if (yPre <= 0.) - continue; - if (TMath::Abs(y) > fConfigRapidityCut) - continue; - } else { - if (yPre == 0.) + /* + // Not sure about this cut. From GammaConv group. Harmless a priori. + if (!(fabs(mctrack.GetEnergy() - mctrack.Pz()) > 0.)) continue; - } - treeWords.fdectyp = mother.getLastDaughterTrackId() - mother.getFirstDaughterTrackId() + 1; // fdectyp: decay type (based on number of daughters). - if (treeWords.fdectyp > 4) - continue; // exclude five or more particles decay - - if (trackID == mctracks.size()) - continue; // no particle left in the list - auto mctrack2 = mctracks[trackID + 1]; - if (!(mctrack2.getMotherTrackId() == mctrack.getMotherTrackId())) - continue; // no matching second electron - if (!(mctrack.getSecondMotherTrackId() == -1)) - continue; // second daughter has more than one mother - if (!(abs(mctrack2.GetPdgCode()) == 11)) - continue; // not an electron - - skipNext = true; // is matching electron --> next particle in list will be skipped - - PxPyPzEVector dau1, dau2, ee; - dau1.SetPxPyPzE(mctrack.Px(), mctrack.Py(), mctrack.Pz(), mctrack.GetEnergy()); - dau2.SetPxPyPzE(mctrack2.Px(), mctrack2.Py(), mctrack2.Pz(), mctrack2.GetEnergy()); - - // create dielectron before resolution effects: - ee = dau1 + dau2; - - // get info of the other particles in the decay: - treeWords.fdau3pdg = 0; - for (Int_t jj = mother.getFirstDaughterTrackId(); jj <= mother.getLastDaughterTrackId(); jj++) { - if (jj == trackID || jj == trackID + 1) { - continue; // first or second electron + // ???? this applied only to first daughter! + Double_t yPre = (mctrack.GetEnergy() + mctrack.Pz()) / (mctrack.GetEnergy() - mctrack.Pz()); + Double_t y = 0.5 * TMath::Log(yPre); + if (fConfigDoRapidityCut) { // Apply rapidity cut on mother consistent with GammaConv group. (??? but it is not applied on mother?) + if (yPre <= 0.) + continue; + if (TMath::Abs(y) > fConfigRapidityCut) + continue; + } else { + if (yPre == 0.) + continue; + }*/ + + treeWords.fdectyp = mother.getLastDaughterTrackId() - mother.getFirstDaughterTrackId() + 1; // fdectyp: decay type (based on number of daughters). + if (treeWords.fdectyp > 4) + continue; // exclude five or more particles decay + + if (trackID == mctracks.size()) + continue; // no particle left in the list + auto mctrack2 = mctracks[trackID + 1]; + if (!(mctrack2.getMotherTrackId() == mctrack.getMotherTrackId())) + continue; // no matching second electron + if (!(mctrack.getSecondMotherTrackId() == -1)) + continue; // second daughter has more than one mother + if (!(abs(mctrack2.GetPdgCode()) == 11)) + continue; // not an electron + + skipNext = true; // is matching electron --> next particle in list will be skipped + + PxPyPzEVector dau1, dau2, ee; + dau1.SetPxPyPzE(mctrack.Px(), mctrack.Py(), mctrack.Pz(), mctrack.GetEnergy()); + dau2.SetPxPyPzE(mctrack2.Px(), mctrack2.Py(), mctrack2.Pz(), mctrack2.GetEnergy()); + + // create dielectron before resolution effects: + ee = dau1 + dau2; + + // get info of the other particles in the decay: + treeWords.fdau3pdg = 0; + for (Int_t jj = mother.getFirstDaughterTrackId(); jj <= mother.getLastDaughterTrackId(); jj++) { + if (jj == trackID || jj == trackID + 1) { + continue; // first or second electron + } + auto mctrack3 = mctracks[jj]; + treeWords.fdau3pdg = abs(mctrack3.GetPdgCode()); } - auto mctrack3 = mctracks[jj]; - treeWords.fdau3pdg = abs(mctrack3.GetPdgCode()); - } - - // get index for histograms - Int_t hindex[3]; - for (Int_t jj = 0; jj < 3; jj++) { - hindex[jj] = -1; - } - switch (mother.GetPdgCode()) { - case 111: - hindex[0] = 0; - break; - case 221: - hindex[0] = 1; - break; - case 331: - hindex[0] = 2; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 22) - hindex[1] = 3; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 223) - hindex[1] = 4; - break; - case 113: - hindex[0] = 5; - break; - case 223: - hindex[0] = 6; - if (treeWords.fdectyp == 2) - hindex[1] = 7; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 111) - hindex[1] = 8; - break; - case 333: - hindex[0] = 9; - if (treeWords.fdectyp == 2) - hindex[1] = 10; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 221) - hindex[1] = 11; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 111) - hindex[1] = 12; - break; - case 443: - hindex[0] = 13; - if (treeWords.fdectyp == 2) - hindex[1] = 14; - if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 22) - hindex[1] = 15; - break; - } - hindex[2] = nInputParticles; - - if (hindex[0] < 0) { - LOGP(error, "hindex[0]<0"); - continue; - } - - // Fill tree words before resolution/acceptance - treeWords.fd1origpt = dau1.Pt(); - treeWords.fd1origp = dau1.P(); - treeWords.fd1origeta = dau1.Eta(); - treeWords.fd1origphi = dau1.Phi(); - treeWords.fd2origpt = dau2.Pt(); - treeWords.fd2origp = dau2.P(); - treeWords.fd2origeta = dau2.Eta(); - treeWords.fd2origphi = dau2.Phi(); - treeWords.feeorigpt = ee.Pt(); - treeWords.feeorigp = ee.P(); - treeWords.feeorigm = ee.M(); - treeWords.feeorigeta = ee.Eta(); - treeWords.feeorigrap = ee.Rapidity(); - treeWords.feeorigphi = ee.Phi(); - if (mctrack.GetPdgCode() > 0) { - treeWords.feeorigphiv = PhiV(dau1, dau2); - } else { - treeWords.feeorigphiv = PhiV(dau2, dau1); - } - - // get the efficiency weight - Int_t effbin = fhwEffpT->FindBin(treeWords.fd1origpt); - treeWords.fwEffpT = fhwEffpT->GetBinContent(effbin); - effbin = fhwEffpT->FindBin(treeWords.fd2origpt); - treeWords.fwEffpT = treeWords.fwEffpT * fhwEffpT->GetBinContent(effbin); - - // Resolution and acceptance - //------------------------- - int ch1 = 1; - int ch2 = 1; - if (mctrack.GetPdgCode() > 0) { - ch1 = -1; - } - if (mctrack2.GetPdgCode() > 0) { - ch2 = -1; - } - dau1 = applySmearingPxPyPzE(ch1, dau1); - dau2 = applySmearingPxPyPzE(ch2, dau2); - - treeWords.fpass = true; - if (dau1.Pt() < fConfigMinPt || dau2.Pt() < fConfigMinPt) - treeWords.fpass = false; // leg pT cut - if (dau1.Pt() > fConfigMaxPt || dau2.Pt() > fConfigMaxPt) - treeWords.fpass = false; // leg pT cut - if (dau1.Vect().Unit().Dot(dau2.Vect().Unit()) > TMath::Cos(fConfigMinOpAng)) - treeWords.fpass = false; // opening angle cut - if (TMath::Abs(dau1.Eta()) > fConfigMaxEta || TMath::Abs(dau2.Eta()) > fConfigMaxEta) - treeWords.fpass = false; - - // get the pair DCA (based in smeared pT) - for (int jj = 0; jj < nbDCAtemplate; jj++) { // loop over DCA templates - if (dau1.Pt() >= DCATemplateEdges[jj] && dau1.Pt() < DCATemplateEdges[jj + 1]) { - treeWords.fd1DCA = fh_DCAtemplates[jj]->GetRandom(); + // get index for histograms + Int_t hindex[3]; + for (Int_t jj = 0; jj < 3; jj++) { + hindex[jj] = -1; } - if (dau2.Pt() >= DCATemplateEdges[jj] && dau2.Pt() < DCATemplateEdges[jj + 1]) { - treeWords.fd2DCA = fh_DCAtemplates[jj]->GetRandom(); + switch (mother.GetPdgCode()) { + case 111: + hindex[0] = 0; + break; + case 221: + hindex[0] = 1; + break; + case 331: + hindex[0] = 2; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 22) + hindex[1] = 3; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 223) + hindex[1] = 4; + break; + case 113: + hindex[0] = 5; + break; + case 223: + hindex[0] = 6; + if (treeWords.fdectyp == 2) + hindex[1] = 7; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 111) + hindex[1] = 8; + break; + case 333: + hindex[0] = 9; + if (treeWords.fdectyp == 2) + hindex[1] = 10; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 221) + hindex[1] = 11; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 111) + hindex[1] = 12; + break; + case 443: + hindex[0] = 13; + if (treeWords.fdectyp == 2) + hindex[1] = 14; + if (treeWords.fdectyp == 3 && treeWords.fdau3pdg == 22) + hindex[1] = 15; + break; } - } - treeWords.fpairDCA = sqrt((pow(treeWords.fd1DCA, 2) + pow(treeWords.fd2DCA, 2)) / 2); - - // Fill tree words after resolution/acceptance - ee = dau1 + dau2; - treeWords.fd1pt = dau1.Pt(); - treeWords.fd1p = dau1.P(); - treeWords.fd1eta = dau1.Eta(); - treeWords.fd1phi = dau1.Phi(); - treeWords.fd2pt = dau2.Pt(); - treeWords.fd2p = dau2.P(); - treeWords.fd2eta = dau2.Eta(); - treeWords.fd2phi = dau2.Phi(); - treeWords.feept = ee.Pt(); - treeWords.feemt = ee.Mt(); - treeWords.feep = ee.P(); - treeWords.feem = ee.M(); - treeWords.feeeta = ee.Eta(); - treeWords.feerap = ee.Rapidity(); - treeWords.feephi = ee.Phi(); - if (mctrack.GetPdgCode() > 0) { - treeWords.feephiv = PhiV(dau1, dau2); - } else { - treeWords.feephiv = PhiV(dau2, dau1); - } - treeWords.fmotherpt = mother.GetPt(); - treeWords.fmotherm = sqrt(pow(mother.GetEnergy(), 2) + pow(mother.GetP(), 2)); - treeWords.fmothermt = sqrt(pow(treeWords.fmotherm, 2) + pow(treeWords.fmotherpt, 2)); - treeWords.fmotherp = mother.GetP(); - treeWords.fmothereta = mother.GetEta(); - treeWords.fmotherphi = mother.GetPhi(); - treeWords.fID = mother.GetPdgCode(); - treeWords.fweight = mctrack.getWeight(); // get particle weight from generator - - // get multiplicity based weight: - int iwbin = fhwMultpT->FindBin(treeWords.fmotherpt); - treeWords.fwMultpT = fhwMultpT->GetBinContent(iwbin); // pT weight - treeWords.fwMultpT2 = fhwMultpT2->GetBinContent(iwbin); // pT weight - double min_mT = fhwMultmT->GetBinLowEdge(1); // consider as minimum valid mT value the edge of the weight histo. - if (treeWords.fmothermt > min_mT) { - iwbin = fhwMultmT->FindBin(treeWords.fmothermt); - treeWords.fwMultmT = fhwMultmT->GetBinContent(iwbin); // mT weight - treeWords.fwMultmT2 = fhwMultmT2->GetBinContent(iwbin); // mT weight - } else { - LOGP(error, "Generated particle with mT < Pion mass cannot be weighted"); - treeWords.fwMultmT = 0.; - treeWords.fwMultmT2 = 0.; - } - // Which ALT weight to use?: - Double_t fwALT = treeWords.fwEffpT; // by default use pt efficiency weight - if (fConfigALTweight == 1) - fwALT = treeWords.fwMultmT; // mT multiplicity weight - if (fConfigALTweight == 11) - fwALT = treeWords.fwMultmT2; // mT multiplicity weight, higher mult - if (fConfigALTweight == 2) - fwALT = treeWords.fwMultpT; // pT multiplicity weight - if (fConfigALTweight == 22) - fwALT = treeWords.fwMultpT2; // pT multiplicity weight, higher mult - - // fill the tree - if (fConfigWriteTTree) { - tree->Fill(); - } + hindex[2] = nInputParticles; - // fill the histograms - if (treeWords.fdectyp < 4) { // why here <4 and before <5 ??? - for (Int_t jj = 0; jj < 3; jj++) { // fill the different hindex -> particles - if (hindex[jj] > -1) { - fmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.fweight); - if (fConfigALTweight == 1 || fConfigALTweight == 11) { - fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmothermt, treeWords.fweight); - } else if (fConfigALTweight == 2 || fConfigALTweight == 22 || fConfigALTweight == 0) { - fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight); - } - fpteevsmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, treeWords.fweight); - fphi_orig[hindex[jj]]->Fill(treeWords.feeorigphi, treeWords.fweight); - frap_orig[hindex[jj]]->Fill(treeWords.feeorigrap, treeWords.fweight); - fmee_orig_wALT[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.fweight * fwALT); - fpteevsmee_orig_wALT[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, treeWords.fweight * fwALT); - if (fConfigALTweight == 1 || fConfigALTweight == 11) { - fmotherpT_orig_wALT[hindex[jj]]->Fill(treeWords.fmothermt, treeWords.fweight * fwALT); - } else if (fConfigALTweight == 2 || fConfigALTweight == 22 || fConfigALTweight == 0) { - fmotherpT_orig_wALT[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight * fwALT); - } - if (treeWords.fpass) { - fmee[hindex[jj]]->Fill(treeWords.feem, treeWords.fweight); - fpteevsmee[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, treeWords.fweight); - fphi[hindex[jj]]->Fill(treeWords.feephi, treeWords.fweight); - frap[hindex[jj]]->Fill(treeWords.feerap, treeWords.fweight); - registry.fill(HIST("DCAeevsmee"), treeWords.feem, treeWords.fpairDCA, treeWords.fweight); - registry.fill(HIST("DCAeevsptee"), treeWords.feept, treeWords.fpairDCA, treeWords.fweight); - fmee_wALT[hindex[jj]]->Fill(treeWords.feem, treeWords.fweight * fwALT); - fpteevsmee_wALT[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, treeWords.fweight * fwALT); - } - } + if (hindex[0] < 0) { + LOGP(error, "hindex[0]<0"); + continue; } - } - if (fConfigDoVirtPh) { - // Virtual photon generation - //------------------------- - // We will generate one virtual photon per histogrammed pion - if (mother.GetPdgCode() == 111) { - // get mass and pt from histos and flat eta and phi - Double_t VPHpT = ffVPHpT->GetRandom(); - Double_t VPHmass = fhKW->GetRandom(); - Double_t VPHeta = -1. + gRandom->Rndm() * 2.; - Double_t VPHphi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm(); - TLorentzVector beam; - beam.SetPtEtaPhiM(VPHpT, VPHeta, VPHphi, VPHmass); - Double_t decaymasses[2] = {(TDatabasePDG::Instance()->GetParticle(11))->Mass(), (TDatabasePDG::Instance()->GetParticle(11))->Mass()}; - TGenPhaseSpace VPHgen; - Bool_t SetDecay; - SetDecay = VPHgen.SetDecay(beam, 2, decaymasses); - if (SetDecay == 0) - LOGP(error, "Decay not permitted by kinematics"); - Double_t VPHweight = VPHgen.Generate(); - // get electrons from the decay - TLorentzVector *decay1, *decay2; - decay1 = VPHgen.GetDecay(0); - decay2 = VPHgen.GetDecay(1); - dau1.SetPxPyPzE(decay1->Px(), decay1->Py(), decay1->Pz(), decay1->E()); - dau2.SetPxPyPzE(decay2->Px(), decay2->Py(), decay2->Pz(), decay2->E()); - - // create dielectron before resolution effects: - ee = dau1 + dau2; - - // get index for histograms - hindex[0] = nInputParticles - 1; - hindex[1] = -1; - hindex[2] = -1; - - // Fill tree words before resolution/acceptance - treeWords.fd1origpt = dau1.Pt(); - treeWords.fd1origp = dau1.P(); - treeWords.fd1origeta = dau1.Eta(); - treeWords.fd1origphi = dau1.Phi(); - treeWords.fd2origpt = dau2.Pt(); - treeWords.fd2origp = dau2.P(); - treeWords.fd2origeta = dau2.Eta(); - treeWords.fd2origphi = dau2.Phi(); - treeWords.feeorigpt = ee.Pt(); - treeWords.feeorigp = ee.P(); - treeWords.feeorigm = ee.M(); - treeWords.feeorigeta = ee.Eta(); - treeWords.feeorigrap = ee.Rapidity(); - treeWords.feeorigphi = ee.Phi(); + // Fill tree words before resolution/acceptance + treeWords.fd1origpt = dau1.Pt(); + treeWords.fd1origp = dau1.P(); + treeWords.fd1origeta = dau1.Eta(); + treeWords.fd1origphi = dau1.Phi(); + treeWords.fd2origpt = dau2.Pt(); + treeWords.fd2origp = dau2.P(); + treeWords.fd2origeta = dau2.Eta(); + treeWords.fd2origphi = dau2.Phi(); + treeWords.feeorigpt = ee.Pt(); + treeWords.feeorigp = ee.P(); + treeWords.feeorigm = ee.M(); + treeWords.feeorigeta = ee.Eta(); + treeWords.feeorigrap = ee.Rapidity(); + treeWords.feeorigphi = ee.Phi(); + if (mctrack.GetPdgCode() > 0) { treeWords.feeorigphiv = PhiV(dau1, dau2); + } else { + treeWords.feeorigphiv = PhiV(dau2, dau1); + } - // get the efficiency weight - Int_t effbin = fhwEffpT->FindBin(treeWords.fd1origpt); - treeWords.fwEffpT = fhwEffpT->GetBinContent(effbin); - effbin = fhwEffpT->FindBin(treeWords.fd2origpt); - treeWords.fwEffpT = treeWords.fwEffpT * fhwEffpT->GetBinContent(effbin); + // get the efficiency weight + Int_t effbin = fhwEffpT->FindBin(treeWords.fd1origpt); + treeWords.fwEffpT = fhwEffpT->GetBinContent(effbin); + effbin = fhwEffpT->FindBin(treeWords.fd2origpt); + treeWords.fwEffpT = treeWords.fwEffpT * fhwEffpT->GetBinContent(effbin); - // Resolution and acceptance - //------------------------- - dau1 = applySmearingPxPyPzE(1, dau1); - dau2 = applySmearingPxPyPzE(-1, dau2); - treeWords.fpass = true; - if (dau1.Pt() < fConfigMinPt || dau2.Pt() < fConfigMinPt) - treeWords.fpass = false; // leg pT cut - if (dau1.Pt() > fConfigMaxPt || dau2.Pt() > fConfigMaxPt) - treeWords.fpass = false; // leg pT cut - if (dau1.Vect().Unit().Dot(dau2.Vect().Unit()) > TMath::Cos(fConfigMinOpAng)) - treeWords.fpass = false; // opening angle cut - if (TMath::Abs(dau1.Eta()) > fConfigMaxEta || TMath::Abs(dau2.Eta()) > fConfigMaxEta) - treeWords.fpass = false; - - treeWords.fpairDCA = 10000.; // ?? - - // Fill tree words after resolution/acceptance - ee = dau1 + dau2; - treeWords.fd1pt = dau1.Pt(); - treeWords.fd1p = dau1.P(); - treeWords.fd1eta = dau1.Eta(); - treeWords.fd1phi = dau1.Phi(); - treeWords.fd2pt = dau2.Pt(); - treeWords.fd2p = dau2.P(); - treeWords.fd2eta = dau2.Eta(); - treeWords.fd2phi = dau2.Phi(); - treeWords.feept = ee.Pt(); - treeWords.feemt = ee.Mt(); - treeWords.feep = ee.P(); - treeWords.feem = ee.M(); - treeWords.feeeta = ee.Eta(); - treeWords.feerap = ee.Rapidity(); - treeWords.feephi = ee.Phi(); - treeWords.feephiv = PhiV(dau1, dau2); - treeWords.fmotherpt = beam.Pt(); - treeWords.fmothermt = sqrt(pow(beam.M(), 2) + pow(beam.Pt(), 2)); - treeWords.fmotherp = beam.P(); - treeWords.fmotherm = beam.M(); - treeWords.fmothereta = beam.Eta(); - treeWords.fmotherphi = beam.Phi(); - treeWords.fID = 0; // set ID to Zero for VPH - treeWords.fweight = VPHweight; - // get multiplicity based weight: - treeWords.fwMultmT = 1; // no weight for photons so far - - // Fill the tree - if (fConfigWriteTTree) { // many parameters not set for photons: d1DCA,fd2DCA, fdectyp,fdau3pdg,fwMultpT,fwMultpT2,fwMultmT2 - tree->Fill(); + // Resolution and acceptance + //------------------------- + int ch1 = 1; + int ch2 = 1; + if (mctrack.GetPdgCode() > 0) { + ch1 = -1; + } + if (mctrack2.GetPdgCode() > 0) { + ch2 = -1; + } + dau1 = applySmearingPxPyPzE(ch1, dau1); + dau2 = applySmearingPxPyPzE(ch2, dau2); + + treeWords.fd1pt = dau1.Pt(); + treeWords.fd1eta = dau1.Eta(); + treeWords.fd2pt = dau2.Pt(); + treeWords.fd2eta = dau2.Eta(); + treeWords.fpass = true; + if (treeWords.fd1pt < fConfigMinPt || treeWords.fd2pt < fConfigMinPt) + treeWords.fpass = false; // leg pT cut + if (treeWords.fd1pt > fConfigMaxPt || treeWords.fd2pt > fConfigMaxPt) + treeWords.fpass = false; // leg pT cut + if (dau1.Vect().Unit().Dot(dau2.Vect().Unit()) > TMath::Cos(fConfigMinOpAng)) + treeWords.fpass = false; // opening angle cut + if (TMath::Abs(treeWords.fd1eta) > fConfigMaxEta || TMath::Abs(treeWords.fd2eta) > fConfigMaxEta) + treeWords.fpass = false; + + // get the pair DCA (based in smeared pT) + for (int jj = 0; jj < nbDCAtemplate; jj++) { // loop over DCA templates + if (dau1.Pt() >= DCATemplateEdges[jj] && dau1.Pt() < DCATemplateEdges[jj + 1]) { + treeWords.fd1DCA = fh_DCAtemplates[jj]->GetRandom(); + } + if (dau2.Pt() >= DCATemplateEdges[jj] && dau2.Pt() < DCATemplateEdges[jj + 1]) { + treeWords.fd2DCA = fh_DCAtemplates[jj]->GetRandom(); } + } + treeWords.fpairDCA = sqrt((pow(treeWords.fd1DCA, 2) + pow(treeWords.fd2DCA, 2)) / 2); + + // Fill tree words after resolution/acceptance + ee = dau1 + dau2; + treeWords.fd1p = dau1.P(); + treeWords.fd1phi = dau1.Phi(); + treeWords.fd2p = dau2.P(); + treeWords.fd2phi = dau2.Phi(); + treeWords.feept = ee.Pt(); + treeWords.feemt = ee.Mt(); + treeWords.feep = ee.P(); + treeWords.feem = ee.M(); + treeWords.feeeta = ee.Eta(); + treeWords.feerap = ee.Rapidity(); + treeWords.feephi = ee.Phi(); + if (mctrack.GetPdgCode() > 0) { + treeWords.feephiv = PhiV(dau1, dau2); + } else { + treeWords.feephiv = PhiV(dau2, dau1); + } + treeWords.fmotherpt = mother.GetPt(); + treeWords.fmotherm = sqrt(pow(mother.GetEnergy(), 2) + pow(mother.GetP(), 2)); + treeWords.fmothermt = sqrt(pow(treeWords.fmotherm, 2) + pow(treeWords.fmotherpt, 2)); + treeWords.fmotherp = mother.GetP(); + treeWords.fmothereta = mother.GetEta(); + treeWords.fmotherphi = mother.GetPhi(); + treeWords.fID = mother.GetPdgCode(); + treeWords.fweight = mctrack.getWeight(); // get particle weight from generator + + // get multiplicity based weight: + int iwbin = fhwMultpT->FindBin(treeWords.fmotherpt); + treeWords.fwMultpT = fhwMultpT->GetBinContent(iwbin); // pT weight + treeWords.fwMultpT2 = fhwMultpT2->GetBinContent(iwbin); // pT weight + double min_mT = fhwMultmT->GetBinLowEdge(1); // consider as minimum valid mT value the edge of the weight histo. + if (treeWords.fmothermt > min_mT) { + iwbin = fhwMultmT->FindBin(treeWords.fmothermt); + treeWords.fwMultmT = fhwMultmT->GetBinContent(iwbin); // mT weight + treeWords.fwMultmT2 = fhwMultmT2->GetBinContent(iwbin); // mT weight + } else { + LOGP(error, "Generated particle with mT < Pion mass cannot be weighted"); + treeWords.fwMultmT = 0.; + treeWords.fwMultmT2 = 0.; + } - // Fill the histograms + // Which ALT weight to use?: + Double_t fwALT = treeWords.fwEffpT; // by default use pt efficiency weight + if (fConfigALTweight == 1) + fwALT = treeWords.fwMultmT; // mT multiplicity weight + if (fConfigALTweight == 11) + fwALT = treeWords.fwMultmT2; // mT multiplicity weight, higher mult + if (fConfigALTweight == 2) + fwALT = treeWords.fwMultpT; // pT multiplicity weight + if (fConfigALTweight == 22) + fwALT = treeWords.fwMultpT2; // pT multiplicity weight, higher mult + + // fill the tree + if (fConfigWriteTTree) { + tree->Fill(); + } + + // fill the histograms + if (treeWords.fdectyp < 4) { // why here <4 and before <5 ??? for (Int_t jj = 0; jj < 3; jj++) { // fill the different hindex -> particles if (hindex[jj] > -1) { - fmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, VPHweight); - fpteevsmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, VPHweight); - fphi_orig[hindex[jj]]->Fill(treeWords.feeorigphi, VPHweight); - frap_orig[hindex[jj]]->Fill(treeWords.feeorigrap, VPHweight); - fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight); + fmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.fweight); + if (fConfigALTweight == 1 || fConfigALTweight == 11) { + fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmothermt, treeWords.fweight); + } else if (fConfigALTweight == 2 || fConfigALTweight == 22 || fConfigALTweight == 0) { + fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight); + } + fpteevsmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, treeWords.fweight); + fphi_orig[hindex[jj]]->Fill(treeWords.feeorigphi, treeWords.fweight); + frap_orig[hindex[jj]]->Fill(treeWords.feeorigrap, treeWords.fweight); + fmee_orig_wALT[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.fweight * fwALT); + fpteevsmee_orig_wALT[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, treeWords.fweight * fwALT); + if (fConfigALTweight == 1 || fConfigALTweight == 11) { + fmotherpT_orig_wALT[hindex[jj]]->Fill(treeWords.fmothermt, treeWords.fweight * fwALT); + } else if (fConfigALTweight == 2 || fConfigALTweight == 22 || fConfigALTweight == 0) { + fmotherpT_orig_wALT[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight * fwALT); + } if (treeWords.fpass) { - fmee[hindex[jj]]->Fill(treeWords.feem, VPHweight); - fpteevsmee[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, VPHweight); - fphi[hindex[jj]]->Fill(treeWords.feephi, VPHweight); - frap[hindex[jj]]->Fill(treeWords.feerap, VPHweight); + fmee[hindex[jj]]->Fill(treeWords.feem, treeWords.fweight); + fpteevsmee[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, treeWords.fweight); + fphi[hindex[jj]]->Fill(treeWords.feephi, treeWords.fweight); + frap[hindex[jj]]->Fill(treeWords.feerap, treeWords.fweight); + registry.fill(HIST("DCAeevsmee"), treeWords.feem, treeWords.fpairDCA, treeWords.fweight); + registry.fill(HIST("DCAeevsptee"), treeWords.feept, treeWords.fpairDCA, treeWords.fweight); + fmee_wALT[hindex[jj]]->Fill(treeWords.feem, treeWords.fweight * fwALT); + fpteevsmee_wALT[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, treeWords.fweight * fwALT); } } } + } - } // mother.pdgCode()==111 - } // fConfigDoVirtPh + if (fConfigDoVirtPh) { + // Virtual photon generation + //------------------------- + // We will generate one virtual photon per histogrammed pion + if (mother.GetPdgCode() == 111) { + // get mass and pt from histos and flat eta and phi + Double_t VPHpT = ffVPHpT->GetRandom(); + Double_t VPHmass = fhKW->GetRandom(); + Double_t VPHeta = -1. + gRandom->Rndm() * 2.; + Double_t VPHphi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm(); + TLorentzVector beam; + beam.SetPtEtaPhiM(VPHpT, VPHeta, VPHphi, VPHmass); + Double_t decaymasses[2] = {(TDatabasePDG::Instance()->GetParticle(11))->Mass(), (TDatabasePDG::Instance()->GetParticle(11))->Mass()}; + TGenPhaseSpace VPHgen; + Bool_t SetDecay; + SetDecay = VPHgen.SetDecay(beam, 2, decaymasses); + if (SetDecay == 0) + LOGP(error, "Decay not permitted by kinematics"); + Double_t VPHweight = VPHgen.Generate(); + // get electrons from the decay + TLorentzVector *decay1, *decay2; + decay1 = VPHgen.GetDecay(0); + decay2 = VPHgen.GetDecay(1); + dau1.SetPxPyPzE(decay1->Px(), decay1->Py(), decay1->Pz(), decay1->E()); + dau2.SetPxPyPzE(decay2->Px(), decay2->Py(), decay2->Pz(), decay2->E()); + + // create dielectron before resolution effects: + ee = dau1 + dau2; + + // get index for histograms + hindex[0] = nInputParticles - 1; + hindex[1] = -1; + hindex[2] = -1; + + // Fill tree words before resolution/acceptance + treeWords.fd1origpt = dau1.Pt(); + treeWords.fd1origp = dau1.P(); + treeWords.fd1origeta = dau1.Eta(); + treeWords.fd1origphi = dau1.Phi(); + treeWords.fd2origpt = dau2.Pt(); + treeWords.fd2origp = dau2.P(); + treeWords.fd2origeta = dau2.Eta(); + treeWords.fd2origphi = dau2.Phi(); + treeWords.feeorigpt = ee.Pt(); + treeWords.feeorigp = ee.P(); + treeWords.feeorigm = ee.M(); + treeWords.feeorigeta = ee.Eta(); + treeWords.feeorigrap = ee.Rapidity(); + treeWords.feeorigphi = ee.Phi(); + treeWords.feeorigphiv = PhiV(dau1, dau2); + + // get the efficiency weight + Int_t effbin = fhwEffpT->FindBin(treeWords.fd1origpt); + treeWords.fwEffpT = fhwEffpT->GetBinContent(effbin); + effbin = fhwEffpT->FindBin(treeWords.fd2origpt); + treeWords.fwEffpT = treeWords.fwEffpT * fhwEffpT->GetBinContent(effbin); + + // Resolution and acceptance + //------------------------- + dau1 = applySmearingPxPyPzE(1, dau1); + dau2 = applySmearingPxPyPzE(-1, dau2); + treeWords.fpass = true; + if (dau1.Pt() < fConfigMinPt || dau2.Pt() < fConfigMinPt) + treeWords.fpass = false; // leg pT cut + if (dau1.Pt() > fConfigMaxPt || dau2.Pt() > fConfigMaxPt) + treeWords.fpass = false; // leg pT cut + if (dau1.Vect().Unit().Dot(dau2.Vect().Unit()) > TMath::Cos(fConfigMinOpAng)) + treeWords.fpass = false; // opening angle cut + if (TMath::Abs(dau1.Eta()) > fConfigMaxEta || TMath::Abs(dau2.Eta()) > fConfigMaxEta) + treeWords.fpass = false; + + treeWords.fpairDCA = 10000.; // ?? + + // Fill tree words after resolution/acceptance + ee = dau1 + dau2; + treeWords.fd1pt = dau1.Pt(); + treeWords.fd1p = dau1.P(); + treeWords.fd1eta = dau1.Eta(); + treeWords.fd1phi = dau1.Phi(); + treeWords.fd2pt = dau2.Pt(); + treeWords.fd2p = dau2.P(); + treeWords.fd2eta = dau2.Eta(); + treeWords.fd2phi = dau2.Phi(); + treeWords.feept = ee.Pt(); + treeWords.feemt = ee.Mt(); + treeWords.feep = ee.P(); + treeWords.feem = ee.M(); + treeWords.feeeta = ee.Eta(); + treeWords.feerap = ee.Rapidity(); + treeWords.feephi = ee.Phi(); + treeWords.feephiv = PhiV(dau1, dau2); + treeWords.fmotherpt = beam.Pt(); + treeWords.fmothermt = sqrt(pow(beam.M(), 2) + pow(beam.Pt(), 2)); + treeWords.fmotherp = beam.P(); + treeWords.fmotherm = beam.M(); + treeWords.fmothereta = beam.Eta(); + treeWords.fmotherphi = beam.Phi(); + treeWords.fID = 0; // set ID to Zero for VPH + treeWords.fweight = VPHweight; + // get multiplicity based weight: + treeWords.fwMultmT = 1; // no weight for photons so far + + // Fill the tree + if (fConfigWriteTTree) { // many parameters not set for photons: d1DCA,fd2DCA, fdectyp,fdau3pdg,fwMultpT,fwMultpT2,fwMultmT2 + tree->Fill(); + } - } // abs(pdgCode())==11 + // Fill the histograms + for (Int_t jj = 0; jj < 3; jj++) { // fill the different hindex -> particles + if (hindex[jj] > -1) { + fmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, VPHweight); + fpteevsmee_orig[hindex[jj]]->Fill(treeWords.feeorigm, treeWords.feept, VPHweight); + fphi_orig[hindex[jj]]->Fill(treeWords.feeorigphi, VPHweight); + frap_orig[hindex[jj]]->Fill(treeWords.feeorigrap, VPHweight); + fmotherpT_orig[hindex[jj]]->Fill(treeWords.fmotherpt, treeWords.fweight); + if (treeWords.fpass) { + fmee[hindex[jj]]->Fill(treeWords.feem, VPHweight); + fpteevsmee[hindex[jj]]->Fill(treeWords.feem, treeWords.feept, VPHweight); + fphi[hindex[jj]]->Fill(treeWords.feephi, VPHweight); + frap[hindex[jj]]->Fill(treeWords.feerap, VPHweight); + } + } + } + + } // mother.pdgCode()==111 + } // fConfigDoVirtPh - } // loop over mctracks + } // abs(pdgCode())==11 - // Clear buffers - eBuff.clear(); - echBuff.clear(); - eweightBuff.clear(); + } // loop over mctracks + + // Clear buffers + eBuff.clear(); + echBuff.clear(); + eweightBuff.clear(); + } } Double_t PhiV(PxPyPzEVector e1, PxPyPzEVector e2) @@ -736,8 +728,8 @@ struct lmeelfcocktail { AxisSpec ptAxis = {fConfigNBinsPtee, fConfigMinPtee, fConfigMaxPtee, "#it{p}_{T,ee} (GeV/c)"}; AxisSpec mAxis = {fConfigNBinsMee, fConfigMinMee, fConfigMaxMee, "#it{m}_{ee} (GeV/c^{2})"}; - AxisSpec phiAxis = {fConfigNBinsPhi, fConfigMinPhi, fConfigMaxPhi, "#it{phi}_{ee}"}; - AxisSpec rapAxis = {fConfigNBinsRap, fConfigMinRap, fConfigMaxRap, "#it{y}_{ee}"}; + AxisSpec phiAxis = {fConfigNBinsPhi, -TMath::TwoPi() / 2, TMath::TwoPi() / 2, "#it{phi}_{ee}"}; + AxisSpec rapAxis = {fConfigNBinsRap, -fConfigMaxAbsRap, fConfigMaxAbsRap, "#it{y}_{ee}"}; registry.add("NEvents", "NEvents", HistType::kTH1F, {{1, 0, 1}}, false); @@ -751,15 +743,8 @@ struct lmeelfcocktail { registry.add("LSmm_orig", "LSmm_orig", HistType::kTH2F, {mAxis, ptAxis}, true); } - // configurable axes crashed the task. Take them out for the moment - // registry.add("DCAeevsmee", "DCAeevsmee", HistType::kTH2F, {{fConfigMBins, "#it{m}_{ee} (GeV/c^{2})"}, {fConfigDCABins, "DCA_{xy}^{ee} (cm)"}}, true); - // registry.add("DCAeevsptee", "DCAeevsptee", HistType::kTH2F, {{fConfigPtBins, "#it{p}_{T,ee} (GeV/c)"}, {fConfigDCABins, "DCA_{xy}^{ee} (cm)"}}, true); - // replace them with hard coded axes - AxisSpec configPtBins = {{0., 0.5, 1, 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5, 7., 7.5, 8.}, "#it{p}_{T,ee} (GeV/c)"}; - AxisSpec configMBins = {{0., 0.08, 0.14, 0.2, 1.1, 2.7, 2.8, 3.2, 5.0}, "#it{m}_{ee} (GeV/c^{2})"}; - AxisSpec configDcaBins = {{0., 0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 3., 4., 5., 7., 10.}, "DCA_{xy}^{ee} (cm)"}; - registry.add("DCAeevsmee", "DCAeevsmee", HistType::kTH2F, {configMBins, configDcaBins}, true); - registry.add("DCAeevsptee", "DCAeevsptee", HistType::kTH2F, {configPtBins, configDcaBins}, true); + registry.add("DCAeevsmee", "DCAeevsmee", HistType::kTH2F, {{fConfigMBins, "#it{m}_{ee} (GeV/c^{2})"}, {fConfigDCABins, "DCA_{xy}^{ee} (cm)"}}, true); + registry.add("DCAeevsptee", "DCAeevsptee", HistType::kTH2F, {{fConfigPtBins, "#it{p}_{T,ee} (GeV/c)"}, {fConfigDCABins, "DCA_{xy}^{ee} (cm)"}}, true); for (auto& particle : fParticleListNames) { fmee.push_back(registry.add(Form("mee_%s", particle.Data()), Form("mee_%s", particle.Data()), HistType::kTH1F, {mAxis}, true)); @@ -987,12 +972,11 @@ struct lmeelfcocktail { { // Build Kroll-wada for virtual photon mass parametrization: Double_t KWmass = 0.; - // Int_t KWnbins = 10000; + Int_t KWnbins = 10000; Float_t KWmin = 2. * eMass; - // Float_t KWmax = 1.1; - Double_t KWbinwidth = (fConfigKWMax - KWmin) / (Double_t)fConfigKWNBins; - fhKW = new TH1F("fhKW", "fhKW", fConfigKWNBins, KWmin, fConfigKWMax); - for (Int_t ibin = 1; ibin <= fConfigKWNBins; ibin++) { + Double_t KWbinwidth = (fConfigKWMax - KWmin) / (Double_t)KWnbins; + fhKW = new TH1F("fhKW", "fhKW", KWnbins, KWmin, fConfigKWMax); + for (Int_t ibin = 1; ibin <= KWnbins; ibin++) { KWmass = KWmin + (Double_t)(ibin - 1) * KWbinwidth + KWbinwidth / 2.0; fhKW->AddBinContent(ibin, 2. * (1. / 137.03599911) / 3. / 3.14159265359 / KWmass * sqrt(1. - 4. * eMass * eMass / KWmass / KWmass) * (1. + 2. * eMass * eMass / KWmass / KWmass)); }