1+ // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+ // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+ // All rights not expressly granted are reserved.
4+ //
5+ // This software is distributed under the terms of the GNU General Public
6+ // License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+ //
8+ // In applying this license CERN does not waive the privileges and immunities
9+ // granted to it by virtue of its status as an Intergovernmental Organization
10+ // or submit itself to any jurisdiction.
11+ // /
12+ // / \file TPCPIDResponse.h
13+ // / \author Annalena Kalteyer <annalena.sophie.kalteyer@cern.ch>
14+ // / \author Christian Sonnabend
15+ // / \brief
16+
17+ #ifndef O2_PID_TPC_RESPONSE_H_
18+ #define O2_PID_TPC_RESPONSE_H_
19+
20+ #include < array>
21+ #include < vector>
22+ #include < cmath>
23+ #include < iostream>
24+ #include " Framework/Logger.h"
25+ // ROOT includes
26+ #include " TFormula.h"
27+ // O2 includes
28+ #include " ReconstructionDataFormats/PID.h"
29+ #include " DataFormatsTPC/BetheBlochAleph.h"
30+
31+ namespace o2 ::pid::tpc
32+ {
33+
34+ // / \brief Class to handle the TPC PID response
35+
36+ class Response
37+ {
38+
39+ public:
40+ Response () = default ;
41+ ~Response () = default ;
42+
43+ // / Setter and Getter for the private parameters
44+ void SetBetheBlochParams (const std::array<float , 5 >& betheBlochParams) { mBetheBlochParams = betheBlochParams; }
45+ void SetResolutionParams (const std::vector<double >& resolutionParams) { mResolutionParams = resolutionParams; }
46+ void SetMIP (const float mip) { mMIP = mip; }
47+ void SetChargeFactor (const float chargeFactor) { mChargeFactor = chargeFactor; }
48+ void SetMultiplicityNormalization (const float multNormalization) { mMultNormalization = multNormalization; }
49+ void SetResolutionParametrization (TFormula* sigmaParametrization) { mSigmaParametrization .reset (sigmaParametrization); }
50+ void SetUseDefaultResolutionParam (const bool useDefault) { mUseDefaultResolutionParam = useDefault; }
51+
52+ const std::array<float , 5 > GetBetheBlochParams () const { return mBetheBlochParams ; }
53+ const std::vector<double > GetResolutionParams () const { return mResolutionParams ; }
54+ const float GetMIP () const { return mMIP ; }
55+ const float GetChargeFactor () const { return mChargeFactor ; }
56+ const float GetMultiplicityNormalization () const { return mMultNormalization ; }
57+ // std::unique_ptr<TFormula> GetResolutionParametrization() { return mSigmaParametrization; }
58+
59+ // / Gets the expected signal of the track
60+ template <typename TrackType>
61+ float GetExpectedSignal (const TrackType& track, const o2::track::PID ::ID id) const ;
62+ // / Gets the expected resolution of the track
63+ template <typename CollisionType, typename TrackType>
64+ float GetExpectedSigma (const CollisionType& collision, const TrackType& trk, const o2::track::PID ::ID id) const ;
65+ // / Gets the number of sigmas with respect the expected value
66+ template <typename CollisionType, typename TrackType>
67+ float GetNumberOfSigma (const CollisionType& collision, const TrackType& trk, const o2::track::PID ::ID id) const ;
68+ // / Gets the deviation to the expected signal
69+ template <typename TrackType>
70+ float GetSignalDelta (const TrackType& trk, const o2::track::PID ::ID id) const ;
71+ // / Gets relative dEdx resolution contribution due to relative pt resolution
72+ float GetRelativeResolutiondEdx (const float p, const float mass, const float charge, const float resol) const ;
73+
74+ void PrintAll () const ;
75+
76+ private:
77+ std::array<float , 5 > mBetheBlochParams = {0.0320981 , 19.9768 , 2.52666e-16 , 2.72123 , 6.08092 };
78+ std::array<float , 2 > mResolutionParamsDefault = {0.07 , 0.0 };
79+ std::vector<double > mResolutionParams = {5.43799e-7 , 0.053044 , 0.667584 , 0.0142667 , 0.00235175 , 1.22482 , 2.3501e-7 , 0.031585 };
80+ float mMIP = 50 .f;
81+ float mChargeFactor = 2 .3f ;
82+ float mMultNormalization = 11000 .;
83+ bool mUseDefaultResolutionParam = true ;
84+ std::unique_ptr<TFormula> mSigmaParametrization {new TFormula (" fSigmaParametrization" , " sqrt(([0]**2)*x[0]+(([1]**2)*(x[2]*[5])*(x[0]/sqrt(1+x[1]**2))**[2])+x[2]*x[3]**2+([4]*x[4])**2 +((x[5]*[6])**2)+(x[5]*(x[0]/sqrt(1+x[1]**2))*[7])**2)" )};
85+
86+ ClassDefNV (Response, 2 );
87+
88+ }; // class Response
89+
90+ // / Get expected Signal of the measurement
91+ template <typename TrackType>
92+ inline float Response::GetExpectedSignal (const TrackType& track, const o2::track::PID ::ID id) const
93+ {
94+ const float bethe = mMIP * o2::tpc::BetheBlochAleph (track.tpcInnerParam () / o2::track::pid_constants::sMasses [id], mBetheBlochParams [0 ], mBetheBlochParams [1 ], mBetheBlochParams [2 ], mBetheBlochParams [3 ], mBetheBlochParams [4 ]) * std::pow ((float )o2::track::pid_constants::sCharges [id], mChargeFactor );
95+ return bethe >= 0 .f ? bethe : 0 .f ;
96+ }
97+
98+ // / Gets the expected resolution of the measurement
99+ template <typename CollisionType, typename TrackType>
100+ inline float Response::GetExpectedSigma (const CollisionType& collision, const TrackType& track, const o2::track::PID ::ID id) const
101+ {
102+ float resolution = 0 .;
103+ if (mUseDefaultResolutionParam ) {
104+ const float reso = track.tpcSignal () * mResolutionParamsDefault [0 ] * ((float )track.tpcNClsFound () > 0 ? std::sqrt (1 . + mResolutionParamsDefault [1 ] / (float )track.tpcNClsFound ()) : 1 .f );
105+ reso >= 0 .f ? resolution = reso : resolution = 0 .f ;
106+ } else {
107+
108+ const double ncl = 159 . / track.tpcNClsFound (); //
109+ const double p = track.tpcInnerParam ();
110+ const double mass = o2::track::pid_constants::sMasses [id];
111+ const double bg = p / mass;
112+ const double dEdx = o2::tpc::BetheBlochAleph ((float )bg, mBetheBlochParams [0 ], mBetheBlochParams [1 ], mBetheBlochParams [2 ], mBetheBlochParams [3 ], mBetheBlochParams [4 ]) * std::pow ((float )o2::track::pid_constants::sCharges [id], mChargeFactor );
113+ const double relReso = GetRelativeResolutiondEdx (p, mass, o2::track::pid_constants::sCharges [id], mResolutionParams [3 ]);
114+
115+ const std::vector<double > values{1 .f / dEdx, track.tgl (), std::sqrt (ncl), relReso, track.signed1Pt (), collision.multTPC () / mMultNormalization };
116+
117+ const float reso = mSigmaParametrization ->EvalPar (values.data (), mResolutionParams .data ()) * dEdx * mMIP ;
118+ reso >= 0 .f ? resolution = reso : resolution = 0 .f ;
119+ }
120+ return resolution;
121+ }
122+
123+ // / Gets the number of sigma between the actual signal and the expected signal
124+ template <typename CollisionType, typename TrackType>
125+ inline float Response::GetNumberOfSigma (const CollisionType& collision, const TrackType& trk, const o2::track::PID ::ID id) const
126+ {
127+ return ((trk.tpcSignal () - GetExpectedSignal (trk, id)) / GetExpectedSigma (collision, trk, id));
128+ }
129+
130+ // / Gets the deviation between the actual signal and the expected signal
131+ template <typename TrackType>
132+ inline float Response::GetSignalDelta (const TrackType& trk, const o2::track::PID ::ID id) const
133+ {
134+ return (trk.tpcSignal () - GetExpectedSignal (trk, id));
135+ }
136+
137+ // // Gets relative dEdx resolution contribution due relative pt resolution
138+ inline float Response::GetRelativeResolutiondEdx (const float p, const float mass, const float charge, const float resol) const
139+ {
140+ const float bg = p / mass;
141+ const float dEdx = o2::tpc::BetheBlochAleph (bg, mBetheBlochParams [0 ], mBetheBlochParams [1 ], mBetheBlochParams [2 ], mBetheBlochParams [3 ], mBetheBlochParams [4 ]) * std::pow (charge, mChargeFactor );
142+ const float deltaP = resol * std::sqrt (dEdx);
143+ const float bgDelta = p * (1 + deltaP) / mass;
144+ const float dEdx2 = o2::tpc::BetheBlochAleph (bgDelta, mBetheBlochParams [0 ], mBetheBlochParams [1 ], mBetheBlochParams [2 ], mBetheBlochParams [3 ], mBetheBlochParams [4 ]) * std::pow (charge, mChargeFactor );
145+ const float deltaRel = std::abs (dEdx2 - dEdx) / dEdx;
146+ return deltaRel;
147+ }
148+
149+ inline void Response::PrintAll () const
150+ {
151+ LOGP (info, " ==== TPC PID response parameters: ====" );
152+ for (int i = 0 ; i < int (mBetheBlochParams .size ()); i++) {
153+ LOGP (info, " BB param [{}] = {}" , i, mBetheBlochParams [i]);
154+ }
155+ LOGP (info, " use default resolution parametrization = {}" , mUseDefaultResolutionParam );
156+ if (mUseDefaultResolutionParam ) {
157+ LOGP (info, " Default Resolution parametrization: " );
158+ for (int i = 0 ; i < int (mResolutionParamsDefault .size ()); i++) {
159+ LOGP (info, " Resolution param [{}] = {}" , i, mResolutionParamsDefault [i]);
160+ }
161+ } else {
162+ for (int i = 0 ; i < int (mResolutionParams .size ()); i++) {
163+ LOGP (info, " Resolution param [{}] = {}" , i, mResolutionParams [i]);
164+ }
165+ }
166+ LOGP (info, " mMIP = {}" , mMIP );
167+ LOGP (info, " mChargeFactor = {}" , mChargeFactor );
168+ LOGP (info, " mMultNormalization = {}" , mMultNormalization );
169+ }
170+
171+ } // namespace o2::pid::tpc
172+
173+ #endif // O2_FRAMEWORK_PIDRESPONSE_H_
0 commit comments