2525#include <TVector3.h>
2626#include "TGrid.h"
2727#include <random>
28+ #include <memory>
2829
2930#include "CCDB/BasicCCDBManager.h"
3031#include "CCDB/CcdbApi.h"
@@ -290,51 +291,55 @@ struct AntinucleiInJets {
290291 }
291292 }
292293
293- void getPerpendicularAxis(TVector3 p, TVector3& u, double sign)
294+ void getPerpendicularAxis(const TVector3& p, TVector3& u, double sign)
294295 {
295- // initialization
296- double ux(0), uy(0), uz(0);
297-
298- // components of vector p
299296 double px = p.X();
300297 double py = p.Y();
301298 double pz = p.Z();
302299
300+ double px2 = px * px;
301+ double py2 = py * py;
302+ double pz2 = pz * pz;
303+ double pz4 = pz2 * pz2;
304+
305+ // px and py are both zero
306+ if (px == 0 && py == 0) {
307+ u.SetXYZ(0, 0, 0);
308+ return;
309+ }
310+
303311 // protection 1
304312 if (px == 0 && py != 0) {
305- uy = -(pz * pz) / py;
306- ux = sign * std::sqrt(py * py - (pz * pz * pz * pz) / (py * py));
307- uz = pz;
308- u.SetXYZ(ux, uy, uz);
313+ double ux = sign * std::sqrt(py2 - pz4 / py2);
314+ double uy = -pz2 / py;
315+ u.SetXYZ(ux, uy, pz);
309316 return;
310317 }
311318
312319 // protection 2
313320 if (py == 0 && px != 0) {
314- ux = -(pz * pz) / px;
315- uy = sign * std::sqrt(px * px - (pz * pz * pz * pz) / (px * px));
316- uz = pz;
317- u.SetXYZ(ux, uy, uz);
321+ double ux = -pz2 / px;
322+ double uy = sign * std::sqrt(px2 - pz4 / px2);
323+ u.SetXYZ(ux, uy, pz);
318324 return;
319325 }
320326
321- // equation parameters
322- double a = px * px + py * py;
323- double b = 2.0 * px * pz * pz;
324- double c = pz * pz * pz * pz - py * py * py * py - px * px * py * py;
327+ // General case
328+ double a = px2 + py2;
329+ double b = 2.0 * px * pz2;
330+ double c = pz4 - py2 * py2 - px2 * py2;
331+
325332 double delta = b * b - 4.0 * a * c;
326333
327- // protection agains delta<0
328- if (delta < 0) {
334+ if (delta < 0 || a == 0) {
335+ LOGP(warn, "Invalid input in getPerpendicularAxis: delta = {}, a = {}", delta, a);
336+ u.SetXYZ(0, 0, 0);
329337 return;
330338 }
331339
332- // solutions
333- ux = (-b + sign * std::sqrt(delta)) / (2.0 * a);
334- uy = (-pz * pz - px * ux) / py;
335- uz = pz;
336- u.SetXYZ(ux, uy, uz);
337- return;
340+ double ux = (-b + sign * std::sqrt(delta)) / (2.0 * a);
341+ double uy = (-pz2 - px * ux) / py;
342+ u.SetXYZ(ux, uy, pz);
338343 }
339344
340345 double getDeltaPhi(double a1, double a2)
@@ -448,19 +453,25 @@ struct AntinucleiInJets {
448453
449454 double getCorrectedPt(double ptRec, TH2* responseMatrix)
450455 {
456+ if (!responseMatrix) {
457+ LOGP(error, "Response matrix is null. Returning uncorrected pt.");
458+ return ptRec;
459+ }
451460
452461 int binX = responseMatrix->GetXaxis()->FindBin(ptRec);
453- TH1D* proj = responseMatrix->ProjectionY("proj", binX, binX);
462+ if (binX < 1 || binX > responseMatrix->GetNbinsX()) {
463+ LOGP(error, "Bin index out of range: binX = {}", binX);
464+ return ptRec; // Return uncorrected pt if bin index is invalid
465+ }
466+ std::unique_ptr<TH1D> proj(responseMatrix->ProjectionY("proj", binX, binX));
454467
455468 // add a protection in case the projection is empty
456469 if (proj->GetEntries() == 0) {
457- delete proj;
458470 return ptRec;
459471 }
460472
461473 double deltaPt = proj->GetRandom();
462474 double ptGen = ptRec + deltaPt;
463- delete proj;
464475
465476 return ptGen;
466477 }
@@ -472,11 +483,12 @@ struct AntinucleiInJets {
472483 LOGP(error, "Could not open the file {}", Form("%s", filepath.Data()));
473484 return;
474485 }
475- responseMatrix = static_cast<TH2F*>( l->FindObject(Form("%s", histoNamePtUnfolding.Data() )));
476- if (!responseMatrix ) {
477- LOGP(error, "Could not open histogram {}", Form("%s", histoNamePtUnfolding.Data()));
486+ TObject* obj = l->FindObject(Form("%s", histoNamePtUnfolding.Data()));
487+ if (!obj || !obj->InheritsFrom(TH2F::Class()) ) {
488+ LOGP(error, "Could not find a valid TH2F histogram {}", Form("%s", histoNamePtUnfolding.Data()));
478489 return;
479490 }
491+ responseMatrix = static_cast<TH2F*>(obj);
480492 LOGP(info, "Opened histogram {}", Form("%s", histoNamePtUnfolding.Data()));
481493 }
482494
0 commit comments