StRoot: StFlowAnalysisMaker/myProjectProfile.C Source File

00001 #include "TH1.h"
00002 #include "TProfile.h"
00003 #include "TProfile2D.h"
00004 #include "iostream.h"
00005 
00006 //absoluteBould =1 make, one can use it like 2.6<|eta|<4
00007 //reversedSign =1 make the addition of profile as  ( negProf *-1. + posProf)
00008 
00009 
00010 TH1D* myProjectionX(TProfile2D* prof2d, double lowBound, double upBound, bool absoluteBound=0, bool reversedSign=0 ){
00011 
00012 
00013 
00014 
00015   TString st(prof2d->GetName());
00016   st +="_X";
00017 
00018 
00019   TH1D* projx = prof2d->ProjectionX(st.Data(), prof2d->GetYaxis()->FindBin(lowBound),prof2d->GetYaxis()->FindBin(upBound));
00020   projx->Reset();
00021 
00022   TH1D* projy = prof2d->ProjectionY();
00023 
00024 
00025 
00026 
00027   for (int i = 0; i<projx->GetNbinsX(); i++) {
00028 
00029     double contentSum=0.;
00030     double entriesSum=0.;
00031     double contentSqrSum=0.;
00032 
00033 
00034     for (int j=0; j<projy->GetNbinsX(); j++) {
00035 
00036 
00037       float BinHighEdge =   (projy->GetBinLowEdge(j) + projy->GetBinWidth(j));
00038       float BinLowEdge = projy->GetBinLowEdge(j);
00039 
00040             if (absoluteBound==0){
00041       if ( BinHighEdge < lowBound) continue;
00042       if ( BinLowEdge  > upBound) continue;
00043             } else if (absoluteBound){
00044 
00045        if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue; 
00046        if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue; 
00047        if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue; 
00048        if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge >  upBound) continue; 
00049        if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
00050             }
00051 
00052 
00053 
00054 
00055 
00056       double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
00057       if (entries==0)  continue;
00058       double content = prof2d->GetBinContent(i,j);
00059       if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.;
00060       double error = prof2d->GetBinError(i,j);
00061 
00062       entriesSum += entries;
00063       contentSum += content * entries;
00064 
00065       if (entries > 1) {
00066       contentSqrSum += entries * (entries * error * error + content * content);
00067       } else if (entries ==1) { //special case to calculate the correct error
00068 
00069         contentSqrSum += content*content;
00070       }
00071     }
00072 
00073     if (entriesSum) {
00074       
00075       projx->SetBinContent(i,contentSum/entriesSum);
00076       projx->SetBinError(i,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
00077 
00078     }
00079 
00080   }
00081 
00082   if (projy) delete projy;
00083   return projx;
00084  
00085 }
00086 
00087 TH1D* myProjectionY(TProfile2D* prof2d, double lowBound, double upBound ,bool absoluteBound=0, bool reversedSign=0  ){
00088 
00089   TString st(prof2d->GetName());
00090   st +="_Y";
00091 
00092 
00093 
00094   TH1D* projy = prof2d->ProjectionY(st.Data(), prof2d->GetXaxis()->FindBin(lowBound),prof2d->GetXaxis()->FindBin(upBound));
00095   projy->Reset();
00096 
00097   TH1D* projx = prof2d->ProjectionX();
00098 
00099 
00100  for (int j=0; j<projy->GetNbinsX(); j++) {
00101 
00102     double contentSum=0.;
00103     double entriesSum=0.;
00104     double contentSqrSum=0.;
00105 
00106      for (int i = 0; i<projx->GetNbinsX(); i++) {
00107 
00108       
00109       float BinHighEdge =   (projx->GetBinLowEdge(i) + projx->GetBinWidth(i));
00110       float BinLowEdge = projx->GetBinLowEdge(i);
00111 
00112             if (absoluteBound==0){
00113       if ( BinHighEdge < lowBound) continue;
00114       if ( BinLowEdge  > upBound) continue;
00115             } else if (absoluteBound){
00116 
00117        if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue; 
00118        if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue; 
00119        if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue; 
00120        if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge >  upBound) continue; 
00121        if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
00122             }
00123 
00124 
00125 
00126 
00127       double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
00128       if (entries==0)  continue;
00129       double content = prof2d->GetBinContent(i,j);
00130       if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.; // for directed flow
00131       double error = prof2d->GetBinError(i,j);
00132 
00133       entriesSum += entries;
00134       contentSum += content * entries;
00135 
00136       if (entries > 1) {
00137       contentSqrSum += entries * (entries * error * error + content * content);
00138       } else if (entries ==1) { //special case to calculate the correct error
00139 
00140         contentSqrSum += content*content;
00141       }
00142     }
00143 
00144     if (entriesSum) {
00145       
00146       projy->SetBinContent(j,contentSum/entriesSum);
00147       projy->SetBinError(j,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
00148 
00149     }
00150 
00151   }
00152  if (projx) delete projx;
00153   return projy;
00154 
00155 }