MpdEmcCluster.cxx 2.8 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283
  1. ////////////////////////////////////////////////////////////////
  2. // //
  3. // MpdEmcHitCluster //
  4. // Cluster production for EMC, v01 //
  5. // Author List : Martemianov M., ITEP, 2017 //
  6. // //
  7. ////////////////////////////////////////////////////////////////
  8. #include <iostream>
  9. #include <numeric>
  10. #include "TMath.h"
  11. #include "MpdEmcCluster.h"
  12. #include "MpdEmcClusterCreation.h"
  13. using namespace std;
  14. using namespace TMath;
  15. // ----- Default constructor -------------------------------------------
  16. MpdEmcCluster::MpdEmcCluster() : fFlag(-1), fE(0.), fTime(0.), fPhi(0.), fRho(0.), fZ(0.){}
  17. // ----- Standard constructor ------------------------------------------
  18. MpdEmcCluster::MpdEmcCluster(Float_t energy, Float_t time, TVector3 pos) :
  19. fFlag(1),
  20. fE(energy),
  21. fTime(time),
  22. fPhi(pos.X()),
  23. fRho(pos.Y()),
  24. fZ(pos.Z()){}
  25. // ----- Destructor ----------------------------------------------------
  26. MpdEmcCluster::~MpdEmcCluster() {
  27. }
  28. // ----- Print -----------------------------------------------------------
  29. void MpdEmcCluster::Print(const Option_t* opt) const {
  30. cout << "MpdEmcCluster: " << endl;
  31. cout << "\tDeposited energy: " << fE <<
  32. "\tMean time: " << fTime <<
  33. " Rho cluster: " << fRho << " Phi cluster: " << fPhi << " Z cluster: " << fZ << endl;
  34. // cout << "\tNumber of tracks in module: " << fNumTracks << endl;
  35. // if (fNumTracks == 1) cout << "PDG code: " << fPDG << " Track ID: " << fTrackID << endl;
  36. }
  37. // -------------------------------------------------------------------------
  38. // Compute claster radius
  39. Float_t MpdEmcCluster::ComputeClusterRadius(TVector3 fPos, vector<Float_t> xC, vector<Float_t> yC,
  40. vector<Float_t> zC, vector<Float_t> eC) {
  41. Float_t phi, theta, valRad, fEnergy, rhoCenter, thetaCenter;
  42. fEnergy = accumulate(eC.begin(), eC.end(), 0.0);
  43. rhoCenter = sqrt(fPos.Y()*fPos.Y()+fPos.Z()*fPos.Z());
  44. thetaCenter = ACos(fPos.Z()/rhoCenter);
  45. valRad = 0.0;
  46. for (UInt_t iHit = 0; iHit < eC.size(); iHit++) {
  47. phi = ATan2(yC[iHit],xC[iHit])*RadToDeg();
  48. if (phi < 0) phi += 360;
  49. theta = ACos(zC[iHit]/sqrt(xC[iHit]*xC[iHit]+yC[iHit]*yC[iHit]+zC[iHit]*zC[iHit]));
  50. valRad += (pow(rhoCenter*(theta-thetaCenter),2)+
  51. pow(fPos.Y()*MpdEmcMath::GetPhiDiff(phi,fPos.X())*DegToRad(),2))*eC[iHit];
  52. /*
  53. valRad += (pow(sqrt(xC[iHit]*xC[iHit]+yC[iHit]*yC[iHit]+zC[iHit]*zC[iHit])*(theta-thetaCenter),2)+
  54. pow(sqrt(xC[iHit]*xC[iHit]+yC[iHit]*yC[iHit])*MpdEmcMath::GetPhiDiff(phi,fPos.X())*DegToRad(),2))*eC[iHit];
  55. */
  56. }
  57. return (sqrt(valRad/fEnergy));
  58. }
  59. ClassImp(MpdEmcCluster)