myProjectProfile.C 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156
  1. #include "TH1.h"
  2. #include "TProfile.h"
  3. #include "TProfile2D.h"
  4. #include "iostream.h"
  5. //absoluteBould =1 make, one can use it like 2.6<|eta|<4
  6. //reversedSign =1 make the addition of profile as ( negProf *-1. + posProf)
  7. TH1D* myProjectionX(TProfile2D* prof2d, double lowBound, double upBound, bool absoluteBound=0, bool reversedSign=0 ){
  8. TString st(prof2d->GetName());
  9. st +="_X";
  10. TH1D* projx = prof2d->ProjectionX(st.Data(), prof2d->GetYaxis()->FindBin(lowBound),prof2d->GetYaxis()->FindBin(upBound));
  11. projx->Reset();
  12. TH1D* projy = prof2d->ProjectionY();
  13. for (int i = 0; i<projx->GetNbinsX(); i++) {
  14. double contentSum=0.;
  15. double entriesSum=0.;
  16. double contentSqrSum=0.;
  17. for (int j=0; j<projy->GetNbinsX(); j++) {
  18. float BinHighEdge = (projy->GetBinLowEdge(j) + projy->GetBinWidth(j));
  19. float BinLowEdge = projy->GetBinLowEdge(j);
  20. if (absoluteBound==0){
  21. if ( BinHighEdge < lowBound) continue;
  22. if ( BinLowEdge > upBound) continue;
  23. } else if (absoluteBound){
  24. if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue;
  25. if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue;
  26. if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue;
  27. if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge > upBound) continue;
  28. if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
  29. }
  30. double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
  31. if (entries==0) continue;
  32. double content = prof2d->GetBinContent(i,j);
  33. if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.;
  34. double error = prof2d->GetBinError(i,j);
  35. entriesSum += entries;
  36. contentSum += content * entries;
  37. if (entries > 1) {
  38. contentSqrSum += entries * (entries * error * error + content * content);
  39. } else if (entries ==1) { //special case to calculate the correct error
  40. contentSqrSum += content*content;
  41. }
  42. }
  43. if (entriesSum) {
  44. projx->SetBinContent(i,contentSum/entriesSum);
  45. projx->SetBinError(i,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
  46. }
  47. }
  48. if (projy) delete projy;
  49. return projx;
  50. }
  51. TH1D* myProjectionY(TProfile2D* prof2d, double lowBound, double upBound ,bool absoluteBound=0, bool reversedSign=0 ){
  52. TString st(prof2d->GetName());
  53. st +="_Y";
  54. TH1D* projy = prof2d->ProjectionY(st.Data(), prof2d->GetXaxis()->FindBin(lowBound),prof2d->GetXaxis()->FindBin(upBound));
  55. projy->Reset();
  56. TH1D* projx = prof2d->ProjectionX();
  57. for (int j=0; j<projy->GetNbinsX(); j++) {
  58. double contentSum=0.;
  59. double entriesSum=0.;
  60. double contentSqrSum=0.;
  61. for (int i = 0; i<projx->GetNbinsX(); i++) {
  62. float BinHighEdge = (projx->GetBinLowEdge(i) + projx->GetBinWidth(i));
  63. float BinLowEdge = projx->GetBinLowEdge(i);
  64. if (absoluteBound==0){
  65. if ( BinHighEdge < lowBound) continue;
  66. if ( BinLowEdge > upBound) continue;
  67. } else if (absoluteBound){
  68. if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue;
  69. if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue;
  70. if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue;
  71. if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge > upBound) continue;
  72. if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
  73. }
  74. double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
  75. if (entries==0) continue;
  76. double content = prof2d->GetBinContent(i,j);
  77. if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.; // for directed flow
  78. double error = prof2d->GetBinError(i,j);
  79. entriesSum += entries;
  80. contentSum += content * entries;
  81. if (entries > 1) {
  82. contentSqrSum += entries * (entries * error * error + content * content);
  83. } else if (entries ==1) { //special case to calculate the correct error
  84. contentSqrSum += content*content;
  85. }
  86. }
  87. if (entriesSum) {
  88. projy->SetBinContent(j,contentSum/entriesSum);
  89. projy->SetBinError(j,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
  90. }
  91. }
  92. if (projx) delete projx;
  93. return projy;
  94. }