GRSISort "v4.0.0.5"
An extension of the ROOT analysis Framework
Loading...
Searching...
No Matches
AngularCorrelationHelper.cxx
Go to the documentation of this file.
2
4{
5 // set up the queue for event mixing by filling it with empty events
6 // we chose 10 here as a reasonable compromise between statistics and computing power/memory needed
7 for(int i = 0; i < fNofMixedEvents; ++i) {
8 fGriffinDeque[slot].emplace_back(new TGriffin);
9 fBgoDeque[slot].emplace_back(new TGriffinBgo);
10 }
11
12 // try and get the cycle length if we have a PPG provided
13 // but only if either the upper or the lower limit of the cycle time is set
14 if(fCycleTimeLow >= 0. || fCycleTimeHigh >= 0.) {
15 if(fPpg != nullptr) {
16 // the ODB cycle length is in microseconds!
18 if(slot == 0) {
19 std::stringstream str;
20 str << "Got ODB cycle length " << fCycleLength << " us = " << fCycleLength / 1e6 << " s" << std::endl;
21 std::cout << str.str();
22 }
23 } else if(slot == 0) {
24 std::stringstream str;
25 str << DRED << "No ppg provided, won't be able to apply cycle cut!" << RESET_COLOR << std::endl;
26 std::cout << str.str();
27 }
28 }
29
30 std::string conditions = fAddback ? "using addback" : "without addback";
31 conditions += fSingleCrystal ? " single crystal" : "";
32 conditions += fFolding ? ", folded around 90^{o}" : "";
33 conditions += fGrouping ? ", grouped" : "";
34
35 for(auto angle = fAngles->begin(); angle != fAngles->end(); ++angle) {
36 int i = std::distance(fAngles->begin(), angle);
37 fH2[slot][Form("AngularCorrelation%d", i)] = new TH2D(Form("AngularCorrelation%d", i), Form("%.1f^{o}: Suppressed #gamma-#gamma %s, |#Deltat_{#gamma-#gamma}| < %.1f", *angle, conditions.c_str(), fPrompt), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
38 fH2[slot][Form("AngularCorrelationBG%d", i)] = new TH2D(Form("AngularCorrelationBG%d", i), Form("%.1f^{o}: Suppressed #gamma-#gamma %s, |#Deltat_{#gamma-#gamma}| = %.1f - %.1f", *angle, conditions.c_str(), fTimeRandomLow, fTimeRandomHigh), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
39 fH2[slot][Form("AngularCorrelationMixed%d", i)] = new TH2D(Form("AngularCorrelationMixed%d", i), Form("%.1f^{o}: Event mixed suppressed #gamma-#gamma %s", *angle, conditions.c_str()), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
40 }
41
42 fH2[slot]["AngularCorrelation"] = new TH2D("AngularCorrelation", Form("Suppressed #gamma-#gamma %s, |#Deltat_{#gamma-#gamma}| < %.1f", conditions.c_str(), fPrompt), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
43 fH2[slot]["AngularCorrelationBG"] = new TH2D("AngularCorrelationBG", Form("Suppressed #gamma-#gamma %s, #Deltat_{#gamma-#gamma} = %.1f - %.1f", conditions.c_str(), fTimeRandomLow, fTimeRandomHigh), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
44 fH2[slot]["AngularCorrelationMixed"] = new TH2D("AngularCorrelationMixed", Form("Event mixed suppressed #gamma-#gamma %s", conditions.c_str()), fBins, fMinEnergy, fMaxEnergy, fBins, fMinEnergy, fMaxEnergy);
45
46 // for the first slot we also write the griffin angles
47 if(slot == 0) {
48 fObject[slot]["GriffinAngles"] = fAngles;
49 }
50}
51
52void AngularCorrelationHelper::Exec(unsigned int slot, TGriffin& fGriffin, TGriffinBgo& fGriffinBgo)
53{
54 bool eventInCycleCut = false;
55 for(auto g1 = 0; g1 < (fAddback ? fGriffin.GetSuppressedAddbackMultiplicity(&fGriffinBgo) : fGriffin.GetSuppressedMultiplicity(&fGriffinBgo)); ++g1) {
56 if(fSingleCrystal && fGriffin.GetNSuppressedAddbackFrags(g1) > 1) continue;
57 auto grif1 = (fAddback ? fGriffin.GetSuppressedAddbackHit(g1) : fGriffin.GetSuppressedHit(g1));
58 if(ExcludeDetector(grif1->GetDetector()) || ExcludeCrystal(grif1->GetArrayNumber())) continue;
59 if(!GoodCycleTime(std::fmod(grif1->GetTime() / 1e3, fCycleLength))) continue;
60
61 // we only get here if at least one hit of the event is within the cycle cut
62 eventInCycleCut = true;
63
64 for(auto g2 = g1 + 1; g2 < (fAddback ? fGriffin.GetSuppressedAddbackMultiplicity(&fGriffinBgo) : fGriffin.GetSuppressedMultiplicity(&fGriffinBgo)); ++g2) {
65 if(fSingleCrystal && fGriffin.GetNSuppressedAddbackFrags(g2) > 1) continue;
66 auto grif2 = (fAddback ? fGriffin.GetSuppressedAddbackHit(g2) : fGriffin.GetSuppressedHit(g2));
67 if(ExcludeDetector(grif2->GetDetector()) || ExcludeCrystal(grif2->GetArrayNumber())) continue;
68 if(!GoodCycleTime(std::fmod(grif2->GetTime() / 1e3, fCycleLength))) continue;
69
70 // skip hits in the same detector when using addback, or in the same crystal when not using addback
71 if(grif1->GetDetector() == grif2->GetDetector() && (fAddback || grif1->GetCrystal() == grif2->GetCrystal())) {
72 continue;
73 }
74
75 // calculate the angle
76 double angle = grif1->GetPosition(fGriffinDistance).Angle(grif2->GetPosition(fGriffinDistance)) * 180. / TMath::Pi();
77 if(angle < fAngles->Rounding()) continue;
78 if(fFolding && angle > 90.) angle = 180. - angle;
79
80 // find the index of the angle
81 auto angleIndex = fAngles->Index(angle);
82 if(angleIndex < 0) continue;
83 // check the timing to see if these are coincident or time-random hits
84 double ggTime = TMath::Abs(grif2->GetTime() - grif1->GetTime());
85 if(ggTime <= fPrompt) {
86 fH2[slot].at("AngularCorrelation")->Fill(grif1->GetEnergy(), grif2->GetEnergy());
87 fH2[slot].at("AngularCorrelation")->Fill(grif2->GetEnergy(), grif1->GetEnergy());
88 fH2[slot].at(Form("AngularCorrelation%d", angleIndex))->Fill(grif1->GetEnergy(), grif2->GetEnergy());
89 fH2[slot].at(Form("AngularCorrelation%d", angleIndex))->Fill(grif2->GetEnergy(), grif1->GetEnergy());
90 } else if(fTimeRandomLow <= ggTime && ggTime <= fTimeRandomHigh) {
91 fH2[slot].at("AngularCorrelationBG")->Fill(grif1->GetEnergy(), grif2->GetEnergy());
92 fH2[slot].at("AngularCorrelationBG")->Fill(grif2->GetEnergy(), grif1->GetEnergy());
93 fH2[slot].at(Form("AngularCorrelationBG%d", angleIndex))->Fill(grif1->GetEnergy(), grif2->GetEnergy());
94 fH2[slot].at(Form("AngularCorrelationBG%d", angleIndex))->Fill(grif2->GetEnergy(), grif1->GetEnergy());
95 }
96 }
97
98 // Event mixing: loop over all stored events for this thread/slot and use them as the "second" gamma ray
99 for(auto l = 0; l < fGriffinDeque[slot].size(); ++l) {
100 auto fLastGriffin = fGriffinDeque[slot][l];
101 auto fLastBgo = fBgoDeque[slot][l];
102 for(auto g2 = 0; g2 < (fAddback ? fLastGriffin->GetSuppressedAddbackMultiplicity(fLastBgo) : fLastGriffin->GetSuppressedMultiplicity(fLastBgo)); ++g2) {
103 auto grif2 = (fAddback ? fLastGriffin->GetSuppressedAddbackHit(g2) : fLastGriffin->GetSuppressedHit(g2));
104 if(ExcludeDetector(grif2->GetDetector()) || ExcludeCrystal(grif2->GetArrayNumber())) continue;
105 if(fSingleCrystal && fLastGriffin->GetNSuppressedAddbackFrags(g2) > 1) continue;
106 if(!GoodCycleTime(std::fmod(grif2->GetTime() / 1e3, fCycleLength))) continue;
107
108 // skip hits in the same detector when using addback, or in the same crystal when not using addback
109 if(grif1->GetDetector() == grif2->GetDetector() && (fAddback || grif1->GetCrystal() == grif2->GetCrystal())) continue;
110
111 double angle = grif1->GetPosition(fGriffinDistance).Angle(grif2->GetPosition(fGriffinDistance)) * 180. / TMath::Pi();
112 if(angle < fAngles->Rounding()) continue;
113 if(fFolding && angle > 90.) angle = 180. - angle;
114
115 // find the index of the angle
116 auto angleIndex = fAngles->Index(angle);
117 if(angleIndex < 0) continue;
118 // no point in checking the time here, the two hits are from different events
119 fH2[slot].at("AngularCorrelationMixed")->Fill(grif1->GetEnergy(), grif2->GetEnergy());
120 fH2[slot].at(Form("AngularCorrelationMixed%d", angleIndex))->Fill(grif1->GetEnergy(), grif2->GetEnergy());
121 }
122 }
123 }
124
125 if(eventInCycleCut) {
126 // Add the current event to the queue, delete the oldest event in the queue, and pop the pointer from the queue.
127 fGriffinDeque[slot].emplace_back(new TGriffin(fGriffin));
128 fBgoDeque[slot].emplace_back(new TGriffinBgo(fGriffinBgo));
129
130 delete fGriffinDeque[slot].front();
131 delete fBgoDeque[slot].front();
132
133 fGriffinDeque[slot].pop_front();
134 fBgoDeque[slot].pop_front();
135 }
136}
#define DRED
Definition Globals.h:18
#define RESET_COLOR
Definition Globals.h:5
void Exec(unsigned int slot, TGriffin &fGriffin, TGriffinBgo &fGriffinBgo)
void CreateHistograms(unsigned int slot) override
Virtual helper function that the user uses to create their histograms.
bool GoodCycleTime(double timeInCycle)
bool ExcludeCrystal(int arraynumber)
std::map< unsigned int, std::deque< TGriffin * > > fGriffinDeque
std::map< unsigned int, std::deque< TGriffinBgo * > > fBgoDeque
std::vector< TGRSIMap< std::string, TH2 * > > fH2
Definition TGRSIHelper.h:49
std::vector< TGRSIMap< std::string, TObject * > > fObject
Definition TGRSIHelper.h:54
TPPG * fPpg
Definition TGRSIHelper.h:56
std::set< double >::iterator end() const
int Index(double angle)
std::set< double >::iterator begin() const
UShort_t GetNSuppressedAddbackFrags(const size_t &idx)
Definition TGriffin.cxx:594
TGriffinHit * GetSuppressedAddbackHit(const int &i)
Definition TGriffin.cxx:539
Short_t GetSuppressedMultiplicity(const TBgo *bgo)
Definition TGriffin.cxx:502
Short_t GetSuppressedAddbackMultiplicity(const TBgo *bgo)
Definition TGriffin.cxx:555
TGriffinHit * GetSuppressedHit(const int &i)
!
Definition TGriffin.cxx:486
Long64_t OdbCycleLength() const
Definition TPPG.h:171