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