1717#include < string>
1818#include < vector>
1919#include < cmath>
20+ #include < set>
21+ #include < utility>
2022
2123#include " Framework/runDataProcessing.h"
2224#include " Framework/AnalysisTask.h"
@@ -55,6 +57,7 @@ using namespace o2::framework::expressions;
5557
5658using myGlobTracks = o2::soa::Join<o2::aod::FullTracks, o2::aod::TrackSelection>;
5759using collisionEvSelIt = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>;
60+ using collisionEvSelItMC = o2::aod::McCollisions;
5861using MCClusters = o2::soa::Join<o2::aod::EMCALMCClusters, o2::aod::EMCALClusters>;
5962
6063struct PhotonIsolationQA {
@@ -78,17 +81,21 @@ struct PhotonIsolationQA {
7881 Configurable<bool > ExoticContribution{" ExoticContribution" , false , " Exotic cluster in the data" };
7982 Configurable<float > minDPhi{" minDPhi" , 0.01 , " Minimum dPhi between track and cluster" };
8083 Configurable<float > Track_matching_Radius{" Track_matching_Radius" , 0.05 , " Radius for which a high energetic track is matched to a cluster" };
84+ Configurable<bool > isMC{" isMC" , true , " should be set to true if the data set is monte carlo" };
8185
8286 Filter PosZFilter = nabs(aod::collision::posZ) < maxPosZ;
87+ Filter PosZFilterMC = nabs(aod::mccollision::posZ) < maxPosZ;
8388 Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == mClusterDefinition ) && (o2::aod::emcalcluster::time >= minTime) && (o2::aod::emcalcluster::time <= maxTime) && (o2::aod::emcalcluster::energy > minClusterEnergy) && (o2::aod::emcalcluster::nCells >= minNCells) && (o2::aod::emcalcluster::nlm <= maxNLM) && (o2::aod::emcalcluster::isExotic == ExoticContribution);
8489 Filter emccellfilter = aod::calo::caloType == 1 ;
8590
8691 using selectedCollisions = soa::Filtered<collisionEvSelIt>;
92+ using selectedMcCollisions = soa::Filtered<collisionEvSelItMC>;
8793 using selectedClusters = soa::Filtered<aod::EMCALClusters>;
8894 using selectedMCClusters = soa::Filtered<MCClusters>;
8995
9096 // Preslices
9197 Preslice<aod::Collisions> collisionsPerBC = aod::collision::bcId;
98+ Preslice<aod::McCollisions> McCollisionsPerBC = aod::mccollision::bcId;
9299 Preslice<aod::Tracks> TracksPercollision = aod::track::collisionId;
93100 Preslice<o2::aod::EMCALMatchedTracks> perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId;
94101 Preslice<selectedMCClusters> ClustersPerCol = aod::emcalcluster::collisionId;
@@ -142,26 +149,41 @@ struct PhotonIsolationQA {
142149 Data_Info.add (" hNCells_NLM_Flag" , " Energy vs NLM" , o2HistType::kTH3F , {{NCells_Axis}, {NLM_Axis}, {SM_Flag_Axis}});
143150
144151 MC_Info.add (" hPosZ" , " Z Position of collision" , o2HistType::kTH1F , {PosZ_Axis});
152+ MC_Info.get <TH1>(HIST (" hPosZ" ))->Sumw2 ();
145153 MC_Info.add (" hNumClusters" , " Number of cluster per collision" , o2HistType::kTH1F , {Num_Cluster_Axis});
154+ MC_Info.get <TH1>(HIST (" hNumClusters" ))->Sumw2 ();
146155 MC_Info.add (" hClusterLocation" , " Location of shower in eta phi plane" , o2HistType::kTH2F , {{Eta_Axis}, {Phi_Axis}});
147156 MC_Info.add (" hEnergy_ShowerShapeLong" , " Energy vs Shower shape long axis" , o2HistType::kTH2F , {{Energy_Axis}, {Shower_Shape_Long_Axis}});
157+ MC_Info.get <TH2>(HIST (" hEnergy_ShowerShapeLong" ))->Sumw2 ();
148158 MC_Info.add (" hEnergy_ShowerShapeShort" , " Energy vs Shower shape short axis" , o2HistType::kTH2F , {{Energy_Axis}, {Shower_Shape_Short_Axis}});
159+ MC_Info.get <TH2>(HIST (" hEnergy_ShowerShapeShort" ))->Sumw2 ();
149160 MC_Info.add (" hEnergy_m02_m20" , " Energy cluster Vs m02 vs m20" , o2HistType::kTH3F , {{Energy_Axis}, {Shower_Shape_Long_Axis}, {Shower_Shape_Short_Axis}});
161+ MC_Info.get <TH3>(HIST (" hEnergy_m02_m20" ))->Sumw2 ();
150162 MC_Info.add (" hEnergy_NCells" , " Energy vs Number of cells in cluster" , o2HistType::kTH2F , {{Energy_Axis}, {NCells_Axis}});
163+ MC_Info.get <TH2>(HIST (" hEnergy_NCells" ))->Sumw2 ();
151164
152165 MC_Info.add (" hEvsNumTracks" , " Energy of cluster vs matched tracks" , o2HistType::kTH2F , {{Energy_Axis}, {Num_Track_Axis}});
166+ MC_Info.get <TH2>(HIST (" hEvsNumTracks" ))->Sumw2 ();
153167 MC_Info.add (" hEvsPtIso" , " Pt_Iso" , o2HistType::kTH2F , {{Energy_Axis}, {PtIso_Axis}});
168+ MC_Info.get <TH2>(HIST (" hEvsPtIso" ))->Sumw2 ();
154169 MC_Info.add (" hRho_Perpen_Cone" , " Energy vs Density of perpendicular cone" , o2HistType::kTH2F , {{Energy_Axis}, {Rho_Axis}});
170+ MC_Info.get <TH2>(HIST (" hRho_Perpen_Cone" ))->Sumw2 ();
155171 MC_Info.add (" hShowerShape" , " Shower shape" , o2HistType::kTH2F , {{Shower_Shape_Long_Axis}, {Shower_Shape_Short_Axis}});
156- MC_Info.add (" hSigmaLongvsPtIso" , " Long shower shape vs Pt_Iso" , o2HistType::kTH2F , {{Shower_Shape_Long_Axis}, {PtIso_Axis}});
172+ MC_Info.get <TH2>(HIST (" hShowerShape" ))->Sumw2 ();
173+ MC_Info.add (" hSigmaLongvsPtIso" , " Long shower shape vs Pt_Iso" , o2HistType::kTH3F , {{Shower_Shape_Long_Axis}, {PtIso_Axis}, {Energy_Axis}});
174+ MC_Info.get <TH3>(HIST (" hSigmaLongvsPtIso" ))->Sumw2 ();
157175 MC_Info.add (" hABCDControlRegion" , " Yield Control Regions" , o2HistType::kTH2F , {{ABCD_Axis}, {Energy_Axis}});
176+ MC_Info.get <TH2>(HIST (" hABCDControlRegion" ))->Sumw2 ();
158177 MC_Info.add (" hClusterEnergy_MCParticleEnergy" , " Energy cluster vs energy particle of cluster" , o2HistType::kTH2F , {{Energy_Axis}, {Energy_Axis}});
178+ MC_Info.get <TH2>(HIST (" hClusterEnergy_MCParticleEnergy" ))->Sumw2 ();
159179 MC_Info.add (" hMotherPDG" , " PDG code of candidate photons mother" , o2HistType::kTH1F , {{2000 , -1000.5 , 999.5 }});
160180 MC_Info.add (" hMotherStatusCode" , " Statuscode of candidate photons mother" , o2HistType::kTH1F , {{400 , -200.5 , 199.5 }});
161181 MC_Info.add (" hMotherStatusCodeVsPDG" , " Statuscode of candidate photons mother" , o2HistType::kTH2F , {{Status_Code_Axis}, {PDG_Axis}});
162182 MC_Info.add (" hCollperBC" , " collisions per BC" , o2HistType::kTH1F , {BC_Axis});
163183 MC_Info.add (" hEnergy_NLM_Flag" , " Energy vs NLM" , o2HistType::kTH3F , {{Energy_Axis}, {NLM_Axis}, {SM_Flag_Axis}});
184+ MC_Info.get <TH3>(HIST (" hEnergy_NLM_Flag" ))->Sumw2 ();
164185 MC_Info.add (" hNCells_NLM_Flag" , " Energy vs NLM" , o2HistType::kTH3F , {{NCells_Axis}, {NLM_Axis}, {SM_Flag_Axis}});
186+ MC_Info.get <TH3>(HIST (" hNCells_NLM_Flag" ))->Sumw2 ();
165187 MC_Info.add (" hPromtPhoton" , " Energy vs m02 vs NCells, PtIso" , o2HistType::kTHnSparseF , {{Energy_Axis}, {Shower_Shape_Long_Axis}, {NCells_Axis}, {PtIso_Axis}});
166188
167189 std::vector<std::string> bin_names = {" A" , " B" , " C" , " D" , " True Bckgr A" };
@@ -248,29 +270,52 @@ struct PhotonIsolationQA {
248270 return Pt_Iso;
249271 }
250272
251- void fillclusterhistos (const auto cluster, HistogramRegistry registry)
273+ void fillclusterhistos (const auto cluster, HistogramRegistry registry, double weight = 1.0 )
252274 {
253275 registry.fill (HIST (" hClusterLocation" ), cluster.eta (), cluster.phi ());
254- registry.fill (HIST (" hEnergy_ShowerShapeLong" ), cluster.energy (), cluster.m02 ());
255- registry.fill (HIST (" hEnergy_ShowerShapeShort" ), cluster.energy (), cluster.m20 ());
256- registry.fill (HIST (" hEnergy_NCells" ), cluster.energy (), cluster.nCells ());
257- registry.fill (HIST (" hEnergy_m02_m20" ), cluster.energy (), cluster.m02 (), cluster.m20 ());
258- registry.fill (HIST (" hShowerShape" ), cluster.m02 (), cluster.m20 ());
276+ if (isMC == true ) {
277+ registry.fill (HIST (" hEnergy_ShowerShapeLong" ), cluster.energy (), cluster.m02 (), weight);
278+ registry.fill (HIST (" hEnergy_ShowerShapeShort" ), cluster.energy (), cluster.m20 (), weight);
279+ registry.fill (HIST (" hEnergy_NCells" ), cluster.energy (), cluster.nCells (), weight);
280+ registry.fill (HIST (" hEnergy_m02_m20" ), cluster.energy (), cluster.m02 (), cluster.m20 (), weight);
281+ registry.fill (HIST (" hShowerShape" ), cluster.m02 (), cluster.m20 (), weight);
282+ } else {
283+ registry.fill (HIST (" hEnergy_ShowerShapeLong" ), cluster.energy (), cluster.m02 ());
284+ registry.fill (HIST (" hEnergy_ShowerShapeShort" ), cluster.energy (), cluster.m20 ());
285+ registry.fill (HIST (" hEnergy_NCells" ), cluster.energy (), cluster.nCells ());
286+ registry.fill (HIST (" hEnergy_m02_m20" ), cluster.energy (), cluster.m02 (), cluster.m20 ());
287+ registry.fill (HIST (" hShowerShape" ), cluster.m02 (), cluster.m20 ());
288+ }
259289 }
260290
261- void fillABCDHisto (HistogramRegistry registry, const auto & cluster, double Pt_iso)
291+ void fillABCDHisto (HistogramRegistry registry, const auto & cluster, double Pt_iso, double weight = 1.0 )
262292 {
263- if ((Pt_iso < 1.5 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
264- registry.fill (HIST (" hABCDControlRegion" ), 0.5 , cluster.energy ());
265- }
266- if ((Pt_iso > 4.0 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
267- registry.fill (HIST (" hABCDControlRegion" ), 1.5 , cluster.energy ());
268- }
269- if ((Pt_iso < 1.5 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
270- registry.fill (HIST (" hABCDControlRegion" ), 2.5 , cluster.energy ());
271- }
272- if ((Pt_iso > 4.0 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
273- registry.fill (HIST (" hABCDControlRegion" ), 3.5 , cluster.energy ());
293+ if (isMC == true ) {
294+ if ((Pt_iso < 1.5 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
295+ registry.fill (HIST (" hABCDControlRegion" ), 0.5 , cluster.energy (), weight);
296+ }
297+ if ((Pt_iso > 4.0 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
298+ registry.fill (HIST (" hABCDControlRegion" ), 1.5 , cluster.energy (), weight);
299+ }
300+ if ((Pt_iso < 1.5 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
301+ registry.fill (HIST (" hABCDControlRegion" ), 2.5 , cluster.energy (), weight);
302+ }
303+ if ((Pt_iso > 4.0 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
304+ registry.fill (HIST (" hABCDControlRegion" ), 3.5 , cluster.energy (), weight);
305+ }
306+ } else {
307+ if ((Pt_iso < 1.5 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
308+ registry.fill (HIST (" hABCDControlRegion" ), 0.5 , cluster.energy ());
309+ }
310+ if ((Pt_iso > 4.0 ) && (cluster.m02 () < 0.3 ) && (cluster.m02 () > 0.1 )) {
311+ registry.fill (HIST (" hABCDControlRegion" ), 1.5 , cluster.energy ());
312+ }
313+ if ((Pt_iso < 1.5 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
314+ registry.fill (HIST (" hABCDControlRegion" ), 2.5 , cluster.energy ());
315+ }
316+ if ((Pt_iso > 4.0 ) && (cluster.m02 () < 2.0 ) && (cluster.m02 () > 0.4 )) {
317+ registry.fill (HIST (" hABCDControlRegion" ), 3.5 , cluster.energy ());
318+ }
274319 }
275320 }
276321
@@ -351,48 +396,48 @@ struct PhotonIsolationQA {
351396 }
352397
353398 // process monte carlo data
354- void processMC (aod::BCs const & bcs, selectedCollisions const & Collisions, selectedMCClusters const & mcclusters, aod::McParticles const &, myGlobTracks const & tracks, o2::aod::EMCALMatchedTracks const & matchedtracks, aod::Calos const &, aod::EMCALClusterCells const & ClusterCells)
399+ void processMC (aod::BCs const & bcs, selectedMcCollisions const & Collisions, selectedMCClusters const & mcclusters, aod::McParticles const &, myGlobTracks const & tracks, o2::aod::EMCALMatchedTracks const & matchedtracks, aod::Calos const &, aod::EMCALClusterCells const & ClusterCells)
355400 {
356401 for (auto bc : bcs) {
357- auto collisionsInBC = Collisions.sliceBy (collisionsPerBC , bc.globalIndex ());
402+ auto collisionsInBC = Collisions.sliceBy (McCollisionsPerBC , bc.globalIndex ());
358403 MC_Info.fill (HIST (" hCollperBC" ), collisionsInBC.size ());
359404 if (collisionsInBC.size () == 1 ) {
360405 for (const auto & Collision : collisionsInBC) {
361- MC_Info.fill (HIST (" hPosZ" ), Collision.posZ ());
406+ MC_Info.fill (HIST (" hPosZ" ), Collision.posZ (), Collision. weight () );
362407 auto ClustersInCol = mcclusters.sliceBy (ClustersPerCol, Collision.globalIndex ());
363408 auto tracksInCol = tracks.sliceBy (TracksPercollision, Collision.globalIndex ());
364409
365410 if (ClustersInCol.size () > 0 ) {
366- MC_Info.fill (HIST (" hNumClusters" ), ClustersInCol.size ());
411+ MC_Info.fill (HIST (" hNumClusters" ), ClustersInCol.size (), Collision. weight () );
367412 }
368413
369414 for (auto & mccluster : ClustersInCol) {
370415 auto tracksofcluster = matchedtracks.sliceBy (perClusterMatchedTracks, mccluster.globalIndex ());
371- fillclusterhistos (mccluster, MC_Info);
372- MC_Info.fill (HIST (" hEvsNumTracks" ), mccluster.energy (), tracksofcluster.size ());
416+ fillclusterhistos (mccluster, MC_Info, Collision. weight () );
417+ MC_Info.fill (HIST (" hEvsNumTracks" ), mccluster.energy (), tracksofcluster.size (), Collision. weight () );
373418
374419 auto CellsInCluster = ClusterCells.sliceBy (CellsPerCluster, mccluster.globalIndex ());
375420 auto [NLM, flag] = CalculateNLM (CellsInCluster);
376- MC_Info.fill (HIST (" hEnergy_NLM_Flag" ), mccluster.energy (), NLM, flag);
377- MC_Info.fill (HIST (" hNCells_NLM_Flag" ), mccluster.nCells (), NLM, flag);
421+ MC_Info.fill (HIST (" hEnergy_NLM_Flag" ), mccluster.energy (), NLM, flag, Collision. weight () );
422+ MC_Info.fill (HIST (" hNCells_NLM_Flag" ), mccluster.nCells (), NLM, flag, Collision. weight () );
378423
379424 if (!track_matching (mccluster, tracksofcluster)) { // no track with significant momentum is matched to cluster
380425 if (NLM <= maxNLM) {
381426 double Pt_Cone = sum_Pt_tracks_in_cone (mccluster, tracksInCol);
382427 double Rho_Perpen_Cone = Rho_Perpendicular_Cone (mccluster, tracksInCol);
383428 double Pt_iso = Pt_Iso (Pt_Cone, Rho_Perpen_Cone);
384429
385- MC_Info.fill (HIST (" hEvsPtIso" ), mccluster.energy (), Pt_iso);
386- MC_Info.fill (HIST (" hRho_Perpen_Cone" ), mccluster.energy (), Rho_Perpen_Cone);
387- MC_Info.fill (HIST (" hSigmaLongvsPtIso" ), mccluster.m02 (), Pt_iso);
388- fillABCDHisto (MC_Info, mccluster, Pt_iso);
430+ MC_Info.fill (HIST (" hEvsPtIso" ), mccluster.energy (), Pt_iso, Collision. weight () );
431+ MC_Info.fill (HIST (" hRho_Perpen_Cone" ), mccluster.energy (), Rho_Perpen_Cone, Collision. weight () );
432+ MC_Info.fill (HIST (" hSigmaLongvsPtIso" ), mccluster.m02 (), Pt_iso, mccluster. energy (), Collision. weight () );
433+ fillABCDHisto (MC_Info, mccluster, Pt_iso, Collision. weight () );
389434
390435 // acces mc true info
391436 auto ClusterParticles = mccluster.mcParticle_as <aod::McParticles>();
392437 bool background = true ;
393438 for (auto & clusterparticle : ClusterParticles) {
394439 if (clusterparticle.pdgCode () == 22 ) {
395- MC_Info.fill (HIST (" hClusterEnergy_MCParticleEnergy" ), mccluster.energy (), clusterparticle.e ());
440+ MC_Info.fill (HIST (" hClusterEnergy_MCParticleEnergy" ), mccluster.energy (), clusterparticle.e (), Collision. weight () );
396441 int first_mother_status_code = getOriginalMotherIndex<aod::McParticles>(clusterparticle);
397442 if (abs (first_mother_status_code) == 23 ) {
398443 background = false ;
@@ -402,7 +447,7 @@ struct PhotonIsolationQA {
402447 }
403448 if (background) {
404449 if ((Pt_iso < 1.5 ) && (mccluster.m02 () < 0.3 ) && (mccluster.m02 () > 0.1 )) {
405- MC_Info.fill (HIST (" hABCDControlRegion" ), 4.5 , mccluster.energy ());
450+ MC_Info.fill (HIST (" hABCDControlRegion" ), 4.5 , mccluster.energy (), Collision. weight () );
406451 }
407452 }
408453 }
0 commit comments