@@ -293,6 +293,7 @@ struct PidFlowPtCorr {
293293 Configurable<bool > cfgOutPutMC1D{" cfgOutPutMC1D" , true , " Fill MC graphs, note that if the processMCgen is open,this MUST be open" };
294294 Configurable<bool > cfgAddPidResponseMatrixHistograms{" cfgAddPidResponseMatrixHistograms" , false , " Add PID response matrix histograms; enable together with processPidResponseMatrix" };
295295 Configurable<bool > cfgAddC22DeltaPtHistograms{" cfgAddC22DeltaPtHistograms" , false , " Add histograms for processDataC22DeltaPt; enable together with processDataC22DeltaPt" };
296+ Configurable<bool > cfgAddPidPtCorrelationProfiles{" cfgAddPidPtCorrelationProfiles" , false , " Add FlowContainer profiles for processPidPtCorrelations; enable together with processPidPtCorrelations" };
296297 Configurable<bool > cfgC22DeltaPtUsePure{" cfgC22DeltaPtUsePure" , false , " true: use PID POI-POI Pure profiles; false: use PID POI-ref and ref-ref profiles" };
297298 Configurable<bool > cfgAddMeanPtCentNbsHistograms{" cfgAddMeanPtCentNbsHistograms" , false , " Add meanptCentNbs TProfile3D histograms" };
298299 // These switches suppress result allocation and filling only. Charged reference
@@ -428,6 +429,7 @@ struct PidFlowPtCorr {
428429 OutputObj<FlowContainer> fFCKa {FlowContainer (" FlowContainerKa" )};
429430 OutputObj<FlowContainer> fFCPr {FlowContainer (" FlowContainerPr" )};
430431 OutputObj<FlowContainer> fFCUnidentified {FlowContainer (" FlowContainerUnidentified" )};
432+ OutputObj<FlowContainer> fFCPidPtCorr {FlowContainer (" FlowContainerPidPtCorr" )};
431433 // end val used for bootstrap
432434
433435 // define global variables
@@ -501,6 +503,7 @@ struct PidFlowPtCorr {
501503 funcProcessQA,
502504 funcProcessPidResponseMatrix,
503505 funcProcessDataC22DeltaPt,
506+ funcProcessPidPtCorrelations,
504507 funcNumber
505508 };
506509
@@ -731,6 +734,9 @@ struct PidFlowPtCorr {
731734 registry.addClone (" hEventCount/processData" , " hEventCount/processSim" );
732735 // processPidResponseMatrix
733736 registry.addClone (" hEventCount/processData" , " hEventCount/processPidResponseMatrix" );
737+ if (switchsOpts.cfgAddPidPtCorrelationProfiles .value ) {
738+ registry.addClone (" hEventCount/processData" , " hEventCount/processPidPtCorrelations" );
739+ }
734740
735741 if (switchsOpts.cfgAddC22DeltaPtHistograms .value ) {
736742 registry.addClone (" hEventCount/processData" , " hEventCount/processDataC22DeltaPt" );
@@ -877,6 +883,30 @@ struct PidFlowPtCorr {
877883 }
878884 // end init fFCPID
879885
886+ if (switchsOpts.cfgAddPidPtCorrelationProfiles .value ) {
887+ auto * oba4PidPtCorr = new TObjArray ();
888+ oba4PidPtCorr->Add (new TNamed (" meanPtPi" , " meanPtPi" ));
889+ oba4PidPtCorr->Add (new TNamed (" meanPtKa" , " meanPtKa" ));
890+ oba4PidPtCorr->Add (new TNamed (" meanPtPr" , " meanPtPr" ));
891+ oba4PidPtCorr->Add (new TNamed (" ptProductPiPi" , " ptProductPiPi" ));
892+ oba4PidPtCorr->Add (new TNamed (" ptPiInPiPiPairs" , " ptPiInPiPiPairs" ));
893+ oba4PidPtCorr->Add (new TNamed (" ptProductKaKa" , " ptProductKaKa" ));
894+ oba4PidPtCorr->Add (new TNamed (" ptKaInKaKaPairs" , " ptKaInKaKaPairs" ));
895+ oba4PidPtCorr->Add (new TNamed (" ptProductPrPr" , " ptProductPrPr" ));
896+ oba4PidPtCorr->Add (new TNamed (" ptPrInPrPrPairs" , " ptPrInPrPrPairs" ));
897+ oba4PidPtCorr->Add (new TNamed (" ptProductPiKa" , " ptProductPiKa" ));
898+ oba4PidPtCorr->Add (new TNamed (" ptPiInPiKaPairs" , " ptPiInPiKaPairs" ));
899+ oba4PidPtCorr->Add (new TNamed (" ptKaInPiKaPairs" , " ptKaInPiKaPairs" ));
900+ oba4PidPtCorr->Add (new TNamed (" ptProductPiPr" , " ptProductPiPr" ));
901+ oba4PidPtCorr->Add (new TNamed (" ptPiInPiPrPairs" , " ptPiInPiPrPairs" ));
902+ oba4PidPtCorr->Add (new TNamed (" ptPrInPiPrPairs" , " ptPrInPiPrPairs" ));
903+ oba4PidPtCorr->Add (new TNamed (" ptProductKaPr" , " ptProductKaPr" ));
904+ oba4PidPtCorr->Add (new TNamed (" ptKaInKaPrPairs" , " ptKaInKaPrPairs" ));
905+ oba4PidPtCorr->Add (new TNamed (" ptPrInKaPrPairs" , " ptPrInKaPrPairs" ));
906+ fFCPidPtCorr ->SetName (" FlowContainerPidPtCorr" );
907+ fFCPidPtCorr ->Initialize (oba4PidPtCorr, axisMultiplicity, cfgFlowNbootstrap);
908+ }
909+
880910 if (switchsOpts.cfgAddMeanPtCentNbsHistograms .value ) {
881911 // TProfile3D axes: event mean pT, centrality, and bootstrap subsample.
882912 if (switchsOpts.cfgOutputCharged .value ) {
@@ -932,6 +962,7 @@ struct PidFlowPtCorr {
932962 registry.add (" hProcessQA/centVsMult" , " cent Vs Mult;Centrality T0C;mulplicity global tracks" , {HistType::kTH2D , {axisMultiplicity, cfgaxisNch}});
933963 registry.add (" hProcessQA/IR" , " " , {HistType::kTH1D , {{100 , 0 , 100 }}});
934964 registry.add (" hProcessQA/Occupacy" , " " , {HistType::kTH1D , {{1000 , 0 , 10000 }}});
965+ registry.add (" hProcessQA/occupancyVsCent" , " Track occupancy versus centrality;Centrality (%);Track occupancy in time range" , {HistType::kTH2D , {axisMultiplicity, {1000 , 0 , 10000 }}});
935966 // end evetn QA
936967 }
937968 // end init QA plot for processQA
@@ -2254,6 +2285,9 @@ struct PidFlowPtCorr {
22542285 case MyFunctionName::funcProcessDataC22DeltaPt:
22552286 registry.fill (HIST (" hEventCount/processDataC22DeltaPt" ), position);
22562287 break ;
2288+ case MyFunctionName::funcProcessPidPtCorrelations:
2289+ registry.fill (HIST (" hEventCount/processPidPtCorrelations" ), position);
2290+ break ;
22572291
22582292 default :
22592293 // LOGF(warning, "could not find event count graph");
@@ -3243,6 +3277,132 @@ struct PidFlowPtCorr {
32433277 }
32443278 PROCESS_SWITCH (PidFlowPtCorr, processDataC22DeltaPt, " Fill c22dmeanpt profiles" , false );
32453279
3280+ void fillPidPtSelfMoments (const char * productName,
3281+ const char * singlePtName,
3282+ const double cent,
3283+ const double rndm,
3284+ const double sumWeight,
3285+ const double sumPt,
3286+ const double sumWeight2,
3287+ const double sumPtWeight2,
3288+ const double sumPt2Weight2)
3289+ {
3290+ const double pairWeight = sumWeight * sumWeight - sumWeight2;
3291+ if (pairWeight <= minVal4Float) {
3292+ return ;
3293+ }
3294+ fFCPidPtCorr ->FillProfile (productName, cent,
3295+ (sumPt * sumPt - sumPt2Weight2) / pairWeight,
3296+ pairWeight, rndm);
3297+ fFCPidPtCorr ->FillProfile (singlePtName, cent,
3298+ (sumWeight * sumPt - sumPtWeight2) / pairWeight,
3299+ pairWeight, rndm);
3300+ }
3301+
3302+ void fillPidPtCrossMoments (const char * productName,
3303+ const char * alphaPtName,
3304+ const char * betaPtName,
3305+ const double cent,
3306+ const double rndm,
3307+ const double alphaWeight,
3308+ const double alphaPtSum,
3309+ const double betaWeight,
3310+ const double betaPtSum)
3311+ {
3312+ if (alphaWeight <= 0 . || betaWeight <= 0 .) {
3313+ return ;
3314+ }
3315+ const double pairWeight = alphaWeight * betaWeight;
3316+ fFCPidPtCorr ->FillProfile (productName, cent,
3317+ alphaPtSum * betaPtSum / pairWeight,
3318+ pairWeight, rndm);
3319+ fFCPidPtCorr ->FillProfile (alphaPtName, cent,
3320+ alphaPtSum / alphaWeight,
3321+ pairWeight, rndm);
3322+ fFCPidPtCorr ->FillProfile (betaPtName, cent,
3323+ betaPtSum / betaWeight,
3324+ pairWeight, rndm);
3325+ }
3326+
3327+ void processPidPtCorrelations (AodCollisions::iterator const & collision,
3328+ aod::BCsWithTimestamps const &,
3329+ AodTracks const & tracks)
3330+ {
3331+ if (!switchsOpts.cfgAddPidPtCorrelationProfiles .value ) {
3332+ return ;
3333+ }
3334+ registry.fill (HIST (" hEventCount/processPidPtCorrelations" ), 0.5 );
3335+ if (tracks.size () < 1 || !collision.sel8 ()) {
3336+ return ;
3337+ }
3338+ registry.fill (HIST (" hEventCount/processPidPtCorrelations" ), 1.5 );
3339+
3340+ const auto cent = getCentrality (collision);
3341+ auto bc = collision.bc_as <aod::BCsWithTimestamps>();
3342+ const double interactionRate = getInteractionRate (bc.timestamp (), bc.runNumber ());
3343+ if (!eventSelected (collision, cent, interactionRate, MyFunctionName::funcProcessPidPtCorrelations)) {
3344+ return ;
3345+ }
3346+
3347+ loadCorrections (bc.timestamp ());
3348+ const double rndm = fRndm ->Rndm ();
3349+ std::array<double , 3 > sumWeight{};
3350+ std::array<double , 3 > sumPt{};
3351+ std::array<double , 3 > sumWeight2{};
3352+ std::array<double , 3 > sumPtWeight2{};
3353+ std::array<double , 3 > sumPt2Weight2{};
3354+
3355+ for (const auto & track : tracks) {
3356+ if (!trackSelectedForFlow (track)) {
3357+ continue ;
3358+ }
3359+ const int pid = getPidConfigurable (track);
3360+ if (pid < MyParticleType::kPion || pid > MyParticleType::kProton || !isWithinPOIPtRange (pid, track.pt ())) {
3361+ continue ;
3362+ }
3363+
3364+ float weight = 1 .f ;
3365+ setParticleNUEWeight (weight, track, cent, pid);
3366+ const std::size_t index = static_cast <std::size_t >(pid - MyParticleType::kPion );
3367+ const double particleWeight = weight;
3368+ const double weight2 = particleWeight * particleWeight;
3369+ const double pt = track.pt ();
3370+ sumWeight[index] += particleWeight;
3371+ sumPt[index] += particleWeight * pt;
3372+ sumWeight2[index] += weight2;
3373+ sumPtWeight2[index] += weight2 * pt;
3374+ sumPt2Weight2[index] += weight2 * pt * pt;
3375+ }
3376+
3377+ constexpr std::size_t kPi = 0 ;
3378+ constexpr std::size_t kKa = 1 ;
3379+ constexpr std::size_t kPr = 2 ;
3380+ if (sumWeight[kPi ] > 0 .) {
3381+ fFCPidPtCorr ->FillProfile (" meanPtPi" , cent, sumPt[kPi ] / sumWeight[kPi ], sumWeight[kPi ], rndm);
3382+ }
3383+ if (sumWeight[kKa ] > 0 .) {
3384+ fFCPidPtCorr ->FillProfile (" meanPtKa" , cent, sumPt[kKa ] / sumWeight[kKa ], sumWeight[kKa ], rndm);
3385+ }
3386+ if (sumWeight[kPr ] > 0 .) {
3387+ fFCPidPtCorr ->FillProfile (" meanPtPr" , cent, sumPt[kPr ] / sumWeight[kPr ], sumWeight[kPr ], rndm);
3388+ }
3389+
3390+ fillPidPtSelfMoments (" ptProductPiPi" , " ptPiInPiPiPairs" , cent, rndm,
3391+ sumWeight[kPi ], sumPt[kPi ], sumWeight2[kPi ], sumPtWeight2[kPi ], sumPt2Weight2[kPi ]);
3392+ fillPidPtSelfMoments (" ptProductKaKa" , " ptKaInKaKaPairs" , cent, rndm,
3393+ sumWeight[kKa ], sumPt[kKa ], sumWeight2[kKa ], sumPtWeight2[kKa ], sumPt2Weight2[kKa ]);
3394+ fillPidPtSelfMoments (" ptProductPrPr" , " ptPrInPrPrPairs" , cent, rndm,
3395+ sumWeight[kPr ], sumPt[kPr ], sumWeight2[kPr ], sumPtWeight2[kPr ], sumPt2Weight2[kPr ]);
3396+
3397+ fillPidPtCrossMoments (" ptProductPiKa" , " ptPiInPiKaPairs" , " ptKaInPiKaPairs" , cent, rndm,
3398+ sumWeight[kPi ], sumPt[kPi ], sumWeight[kKa ], sumPt[kKa ]);
3399+ fillPidPtCrossMoments (" ptProductPiPr" , " ptPiInPiPrPairs" , " ptPrInPiPrPairs" , cent, rndm,
3400+ sumWeight[kPi ], sumPt[kPi ], sumWeight[kPr ], sumPt[kPr ]);
3401+ fillPidPtCrossMoments (" ptProductKaPr" , " ptKaInKaPrPairs" , " ptPrInKaPrPairs" , cent, rndm,
3402+ sumWeight[kKa ], sumPt[kKa ], sumWeight[kPr ], sumPt[kPr ]);
3403+ }
3404+ PROCESS_SWITCH (PidFlowPtCorr, processPidPtCorrelations, " Calculate self-contained PID mean-pT correlations" , false );
3405+
32463406 /* *
32473407 * @brief Run the flow calculation on generated MC particles for a closure test.
32483408 * @note Reconstructed collisions are only used to obtain the centrality. The selection is
@@ -3675,6 +3835,7 @@ struct PidFlowPtCorr {
36753835 registry.fill (HIST (" hProcessQA/centVsMult" ), cent, tracks.size ());
36763836 registry.fill (HIST (" hProcessQA/IR" ), interactionRate);
36773837 registry.fill (HIST (" hProcessQA/Occupacy" ), collision.trackOccupancyInTimeRange ());
3838+ registry.fill (HIST (" hProcessQA/occupancyVsCent" ), cent, collision.trackOccupancyInTimeRange ());
36783839 // end event qa
36793840
36803841 // i dont want to fill event count again as it's filled in other 5 function
0 commit comments