PWGLF/totEt/macros/emEt/PlotEmEtVer2.C

   1 Int_t colors[] = {TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
   2                   TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
   3                   TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack,
   4                   TColor::kRed, TColor::kOrange, TColor::kGreen+3, TColor::kBlue, TColor::kBlack};
   5 Int_t markers[] = {20,21,22,23,33, 24,25,26,32,27, 20,21,22,23,33, 24,25,26,32,27};
   6 Int_t CutSet = 0;//Defaults: 250 MeV for PHOS and 300 MeV for EMCal
   7 //1:  350 MeV for both
   8 void SetStyles(TGraph *graph, Int_t marker, Int_t color){
   9   graph->SetMarkerStyle(marker);
  10   graph->SetMarkerColor(color);
  11   graph->SetLineColor(color);
  12   graph->SetMarkerSize(1.5);
  13 }
  14 void WriteLatex();
  15 Float_t finalemetCorrEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  16 Float_t finalemetCorrPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  17 Float_t finalemetErrorEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  18 Float_t finalemetErrorPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  19 Float_t finaltotaletCorrEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  20 Float_t finaltotaletCorrPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  21 Float_t finaltotaletErrorEmcal[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  22 Float_t finaltotaletErrorPhos[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  23 //neutron
  24 Float_t GetMostLikelyValue(TH1 *histo){
  25   Float_t max = 0;
  26   Float_t maxx=0;
  27   for(Int_t bin=1;bin<histo->GetNbinsX();bin++){
  28     if(histo->GetBinContent(bin)>max){
  29       max = histo->GetBinContent(bin);
  30       maxx = histo->GetBinCenter(bin);
  31     }
  32   }
  33   return maxx;
  34 }
  35
  36 Float_t npartShort[10] =    {382.8,329.7,260.5,186.4,128.9, 85,52.8,30.0,15.8,7.48};
  37 Float_t npartErrShort[10] = {    6,    6,  4.4,  3.9,  3.3,2.6,   2, 1.3, 0.6,0.29};
  38 Float_t npart[20] = {382.7, 329.4, 281.2, 239, 202.1, 169.5, 141, 116, 94.11, 75.3, 59.24, 45.58, 34.33, 25.21, 17.96, 12.58, 8.812, 6.158, 4.376, 3.064};
  39 Float_t npartErr[20] = {3, 4.3, 4.1, 3.5, 3.3, 3.3, 3.1, 2.8, 2.6, 2.3, 1.8, 1.4, 1.1, 0.87, 0.66, 0.45, 0.26, 0.19, 0.1, 0.059};
  40
  41 //========================Charged Pion Reference========================================
  42 //Arrays for defining comparison plots
  43 Float_t pionPlusEt[10] = {360.7,298.3,223.8,149.9,96.1, 58.1,32.4,16.4,7.3,2.7};
  44 Float_t pionMinusEt[10] ={363.7,300.4,225.4,150.5,96.6, 58.4,32.5,16.5,7.4,2.8};
  45 Float_t pionEtError[10] = {19.3, 15.3,11.3 ,7.5  , 4.8,  2.9, 1.6, 0.8,0.4,0.1};
  46 Float_t pionEt[10] = {0,0,0,0,0, 0,0,0,0,0};
  47 Float_t ypion[10]  = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  48 Float_t ypionerr[10]  = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  49 Float_t emEtFromPions[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  50 Float_t emEtFromPionsErr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  51 Float_t emEtFromPionsPerNpart[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  52 Float_t emEtFromPionsPerNpartErr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  53 Float_t npartAlt1[20],npartAlt2[20],npartAlt3[20];
  54
  55 TGraphErrors *GetPionEtGraph(){
  56   for(int i=0;i<10;i++){
  57     pionEt[i] = (pionPlusEt[i]+pionMinusEt[i])/2.0;
  58     emEtFromPions[i] = pionEt[i]*1.085;
  59     emEtFromPionsErr[i] = emEtFromPions[i]*TMath::Sqrt(TMath::Power(0.030/1.085,2)+TMath::Power(pionEtError[i]/pionEt[i],2));
  60     ypion[i] = pionEt[i]/(npartShort[i]/2);
  61     ypionerr[i] = pionEtError[i]/(npartShort[i]/2);
  62     emEtFromPionsPerNpart[i] = emEtFromPions[i]/(npartShort[i]/2);
  63     emEtFromPionsPerNpartErr[i] = emEtFromPionsErr[i]/(npartShort[i]/2);
  64     npartAlt1[i] = npartShort[i]+2;
  65     npartAlt2[i] = npartShort[i]-2;
  66     npartAlt3[i] = npartShort[i]+4;
  67   }
  68   TGraphErrors *gr2 = new TGraphErrors(10,npartShort,ypion,npartErrShort,ypionerr);
  69   gr2->GetYaxis()->SetTitle("dE_{T}/d#eta#frac{1}{0.5*N_{part}} [GeV]");
  70   gr2->GetXaxis()->SetTitle("N_{part}");
  71   gr2->SetTitle("");
  72   gr2->GetXaxis()->SetRangeUser(0, 400);
  73   SetStyles(gr2,30,TColor::kBlue);
  74
  75   return gr2;
  76 }
  77
  78 TGraphErrors *GetPionEmEtGraph(){
  79     TGraphErrors *gr3 = new TGraphErrors(10,npartAlt3,emEtFromPionsPerNpart,npartErrShort,emEtFromPionsPerNpartErr);
  80     SetStyles(gr3,29,TColor::kBlue);
  81     return gr3;
  82
  83 }
  84 //========================Reading in corrections========================================
  85 Float_t nonLinError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  86 Float_t nonLinErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  87 //Float_t signalFraction[20] = {0.4,0.4,0.4,0.4,0.4, 0.4,0.4,0.4,0.4,0.4, 0.4,0.4,0.4,0.4,0.4, 0.4,0.4,0.4,0.4,0.4};
  88 Float_t signalFraction[20] = {0.3,0.3,0.3,0.3,0.3, 0.3,0.3,0.3,0.3,0.3, 0.3,0.3,0.3,0.3,0.3, 0.3,0.3,0.3,0.3,0.3};
  89 Float_t signalFractionError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  90 //Float_t averageEfficiency[20] = {0.5,0.5,0.5,0.5,0.5, 0.5,0.5,0.5,0.5,0.5, 0.5,0.5,0.5,0.5,0.5, 0.5,0.5,0.5,0.5,0.5};
  91 Float_t averageEfficiency[20] = {1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0,1.0,1.0};
  92 Float_t averageEfficiencyError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  93 Float_t efficiencyError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  94 Float_t efficiencyErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  95 Float_t neutronCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  96 Float_t neutronError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
  97 Float_t neutronErrorShort[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
  98 Float_t neutronCorrShort[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
  99 Float_t neutronErrorPerNChShort[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 100 Float_t neutronCorrPerNChShort[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 101 Float_t neutronErrorPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 102 Float_t neutronCorrPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 103 Float_t secondaryErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 104 Float_t secondaryCorrShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 105 Float_t secondaryCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 106 Float_t secondaryError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 107 Float_t secondaryCorrPerNChShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 108 Float_t secondaryErrorPerNChShort[10] =  {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 109 Float_t secondaryCorrPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 110 Float_t secondaryErrorPerNCh[20] =  {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 111 Float_t kaonError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 112 Float_t kaonCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 113 Float_t kaonErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 114 Float_t kaonCorrShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 115 Float_t kaonErrorPerNChShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 116 Float_t kaonCorrPerNChShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 117 Float_t kaonErrorPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 118 Float_t kaonCorrPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 119 Float_t hadErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 120 Float_t hadCorrShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 121 Float_t hadError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 122 Float_t hadCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 123 Float_t hadronErrorPerNChShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 124 Float_t hadronCorrPerNChShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 125 Float_t hadronErrorPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 126 Float_t hadronCorrPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 127 Float_t minEtErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 128 Float_t minEtCorrShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 129 Float_t minEtError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 130 Float_t minEtCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 131
 132
 133 Float_t neutronCorrNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 134 Float_t neutronErrorNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 135 Float_t neutronErrorShortNoEffCorr[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 136 Float_t neutronCorrShortNoEffCorr[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 137 Float_t neutronErrorPerNChShortNoEffCorr[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 138 Float_t neutronCorrPerNChShortNoEffCorr[11] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0};
 139 Float_t neutronErrorPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 140 Float_t neutronCorrPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 141 Float_t secondaryErrorShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 142 Float_t secondaryCorrShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 143 Float_t secondaryCorrNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 144 Float_t secondaryErrorNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 145 Float_t secondaryCorrPerNChShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 146 Float_t secondaryErrorPerNChShortNoEffCorr[10] =  {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 147 Float_t secondaryCorrPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 148 Float_t secondaryErrorPerNChNoEffCorr[20] =  {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 149 Float_t kaonErrorNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 150 Float_t kaonCorrNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 151 Float_t kaonErrorShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 152 Float_t kaonCorrShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 153 Float_t kaonErrorPerNChShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 154 Float_t kaonCorrPerNChShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 155 Float_t kaonErrorPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 156 Float_t kaonCorrPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 157 Float_t hadErrorShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 158 Float_t hadCorrShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 159 Float_t hadErrorNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 160 Float_t hadCorrNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 161 Float_t hadronErrorPerNChShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 162 Float_t hadronCorrPerNChShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 163 Float_t hadronErrorPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 164 Float_t hadronCorrPerNChNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 165
 166
 167
 168
 169 Float_t kaonYield[10] = {109,90.5,68,46,30,18.2,10.2,5.1,2.3,0.855};
 170 Float_t kaonYieldStatErr[10] = {0.3,0.2,0.1,0.1,0.1, 0.06,0.04,0.03,0.02,0.01};
 171 Float_t kaonYieldSysErr[10] = {9,7,5,4,2, 1.5,0.8,0.4,0.2,0.09};
 172 Float_t kaonYieldTotErr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 173 Float_t kaonYieldPerNCh[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 174 Float_t kaonYieldPerNChErr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 175 Float_t kaonEtPerNCh[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 176 Float_t kaonEtPerNChErr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 177 Double_t kaonPlusEt[2][10] = {{91.7712,75.8971,56.841,37.6962,23.9923,14.255,7.73469,3.74477,1.60505,0.578278},{6.61811,5.39337,3.9978,2.6337,1.67854,1.01849,0.557879,0.278199,0.125057,0.0592682}};
 178 Double_t kaonMinusEt[2][10] = {{90.4723,74.9444,55.9463,37.286,23.6591,14.0413,7.63067,3.69337,1.59219,0.571019},{7.01588,5.76588,4.20933,2.80388,1.77983,1.06934,0.588003,0.292737,0.138191,0.0600075}};
 179
 180 Float_t averageHadEnergy = -1;
 181 Float_t averageHadEnergyError = -1;
 182
 183 void ReadMinEtCorrections(){
 184   cout<<"Reading in min et corrections..."<<endl;
 185   string inline;
 186   float value = 0;
 187   float error = 0;
 188   float err = 0;
 189   int i=0;
 190   //file names like /tmp/MinEtEmcalShortCut7.dat
 191   //1 = 100 MeV/c
 192   //2 = 150 MeV/c
 193   //3 = 200 MeV/c
 194   //4 = 250 MeV/c
 195   //5 = 300 MeV/c
 196   //6 = 350 MeV/c
 197   //7 = 400 MeV/c
 198   //8 = 450 MeV/c
 199   //9 = 500 MeV/c
 200   //10 = 550 MeV/c
 201   TString cutstring = "";
 202   if(CutSet==1) cutstring = "Cut6";
 203   if(CutSet==2) cutstring = "Cut7";
 204   TString minetInfileNameShort = "MinEt"+detector+"Short"+cutstring+".dat";
 205   cout<<"Reading "<<minetInfileNameShort.Data()<<endl;
 206   ifstream myminetfileShort (minetInfileNameShort.Data());
 207   if (myminetfileShort.is_open()){
 208     while ( myminetfileShort.good() )
 209       {
 210         getline (myminetfileShort,inline);
 211         istringstream tmp(inline);
 212         tmp >> value;
 213         tmp >> error;
 214         if(i<10){
 215           minEtCorrShort[i] = value;
 216           minEtErrorShort[i] = error;
 217         }
 218         //cout<<"min et corr cb "<<i<<" "<< minEtCorrShort[i]<<" +/- "<<minEtErrorShort[i]<<endl;
 219         i++;
 220       }
 221     myminetfileShort.close();
 222   }
 223   //doing the linear interpolation between the two data points we have
 224   TGraphErrors *graphMinEtCorrectionShort = GetMinEtCorrectionGraphShort();
 225   int shortbin = 0;
 226   for(int i=0;i<19;i++){
 227     //cout<<"Long bin "<<i<<" short bin "<<shortbin;
 228     if(i<2){//we have exact numbers so we don't need to interpolate
 229       minEtCorr[i] = minEtCorrShort[i];
 230       minEtError[i] = minEtErrorShort[i];
 231       shortbin++;
 232     }
 233     else{
 234       minEtCorr[i] = graphMinEtCorrectionShort->Eval(trackmultiplicity[i]);
 235       int altbin = shortbin-1;
 236       if(i%2==1){altbin = shortbin+1;}
 237       //cout<<" altbin "<<altbin;
 238       if(minEtErrorShort[shortbin]>minEtErrorShort[altbin]) minEtError[i] = minEtErrorShort[shortbin];
 239       else{minEtError[i] = minEtErrorShort[altbin];}
 240       if(i%2==1 && shortbin<10){shortbin++;}
 241     }
 242     // cout<<"min et corr cb "<<i<<" "<< minEtCorr[i]<<" +/- "<<minEtError[i]<<endl;
 243     //cout<<endl;
 244   }
 245   delete graphMinEtCorrectionShort;
 246
 247 }
 248 void ReadInNeutronCorrections(){
 249   cout<<"Reading in neutron corrections..."<<endl;
 250   TString neutronInfileName = "Neutrons"+detector+".dat";
 251   ifstream myneutronfile (neutronInfileName.Data());
 252   string inline;
 253   float value = 0;
 254   float error = 0;
 255   int i=0;
 256   if (myneutronfile.is_open()){
 257     while ( myneutronfile.good() )
 258       {
 259         getline (myneutronfile,inline);
 260         istringstream tmp(inline);
 261         tmp >> value;
 262         tmp >> error;
 263         if(i<20){
 264           neutronCorr[i] = value;
 265           neutronError[i] = error;
 266           if(trackmultiplicity[i]>0){
 267             neutronCorrPerNCh[i] = value/(trackmultiplicity[i]);
 268             neutronErrorPerNCh[i] = error/(trackmultiplicity[i]);
 269           }
 270         }
 271
 272         // cout<<"neutroncorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 273         i++;
 274       }
 275     myneutronfile.close();
 276
 277   }
 278   TString neutronInfileNameShort = "Neutrons"+detector+"Short.dat";
 279   ifstream myneutronfileShort (neutronInfileNameShort.Data());
 280   i=0;
 281   if (myneutronfileShort.is_open()){
 282     while ( myneutronfileShort.good() )
 283       {
 284         getline (myneutronfileShort,inline);
 285         istringstream tmp(inline);
 286         tmp >> value;
 287         tmp >> error;
 288         if(i<10){
 289           neutronCorrShort[i] = value;
 290           neutronErrorShort[i] = error;
 291           neutronCorrPerNChShort[i] = value/(trackmultiplicityShort[i]);
 292           neutronErrorPerNChShort[i] = error/(trackmultiplicityShort[i]);
 293         }
 294         //cout<<"neutroncorr cb "<<i<<" "<<neutronCorrShort[i]<<" +/- "<<neutronErrorShort[i]<<endl;
 295         i++;
 296       }
 297     myneutronfileShort.close();
 298   }
 299   //Begin reading in no eff corr corrections
 300   TString neutronInfileName3 = "Neutrons"+detector+"NoEffCorr.dat";
 301   ifstream myneutronfile3 (neutronInfileName3.Data());
 302   value = 0;
 303   error = 0;
 304   i=0;
 305   if (myneutronfile3.is_open()){
 306     while ( myneutronfile3.good() )
 307       {
 308         getline (myneutronfile3,inline);
 309         istringstream tmp(inline);
 310         tmp >> value;
 311         tmp >> error;
 312         if(i<20){
 313           neutronCorrNoEffCorr[i] = value;
 314           neutronErrorNoEffCorr[i] = error;
 315           if(trackmultiplicity[i]>0){
 316             neutronCorrPerNChNoEffCorr[i] = value/(trackmultiplicity[i]);
 317             neutronErrorPerNChNoEffCorr[i] = error/(trackmultiplicity[i]);
 318           }
 319         }
 320
 321         //cout<<"neutroncorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 322         i++;
 323       }
 324     myneutronfile3.close();
 325
 326   }
 327   TString neutronInfileNameShort3 = "Neutrons"+detector+"NoEffCorrShort.dat";
 328   ifstream myneutronfile3Short (neutronInfileNameShort3.Data());
 329   i=0;
 330   if (myneutronfile3Short.is_open()){
 331     while ( myneutronfile3Short.good() )
 332       {
 333         getline (myneutronfile3Short,inline);
 334         istringstream tmp(inline);
 335         tmp >> value;
 336         tmp >> error;
 337         if(i<10){
 338           neutronCorrShortNoEffCorr[i] = value;
 339           neutronErrorShortNoEffCorr[i] = error;
 340           neutronCorrPerNChShortNoEffCorr[i] = value/(trackmultiplicityShort[i]);
 341           neutronErrorPerNChShortNoEffCorr[i] = error/(trackmultiplicityShort[i]);
 342         }
 343         //cout<<"neutroncorr cb "<<i<<" "<<neutronCorrShortNoEffCorr[i]<<" +/- "<<neutronErrorShortNoEffCorr[i]<<endl;
 344         i++;
 345       }
 346     myneutronfile3Short.close();
 347   }
 348
 349
 350 }
 351 void ReadInSecondaryCorrections(){
 352   cout<<"Reading in secondary corrections..."<<endl;
 353
 354     TString secondaryInfileName = "Secondaries"+detector+".dat";
 355     ifstream mysecondaryfile (secondaryInfileName.Data());
 356     string inline;
 357     float value = 0;
 358     float error = 0;
 359     int i=0;
 360     if (mysecondaryfile.is_open()){
 361       while ( mysecondaryfile.good() )
 362         {
 363           getline (mysecondaryfile,inline);
 364           istringstream tmp(inline);
 365           tmp >> value;
 366           tmp >> error;
 367           if(i<20){
 368             secondaryCorr[i] = value;
 369             secondaryError[i] = error;
 370             if(trackmultiplicity[i]>0){
 371               secondaryCorrPerNCh[i] = value/(trackmultiplicity[i])/5.0;
 372               secondaryErrorPerNCh[i] = error/(trackmultiplicity[i])/5.0;
 373             }
 374           }
 375           //cout<<"secondarycorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 376           i++;
 377         }
 378         mysecondaryfile.close();
 379     }
 380     TString secondaryShortInfileName = "Secondaries"+detector+"Short.dat";
 381     ifstream mysecondaryShortfile (secondaryShortInfileName.Data());
 382     i=0;
 383     if (mysecondaryShortfile.is_open()){
 384       while ( mysecondaryShortfile.good() && i<10 )
 385         {
 386           getline (mysecondaryShortfile,inline);
 387           istringstream tmp(inline);
 388           tmp >> value;
 389           tmp >> error;
 390           if(i<10){
 391             secondaryCorrShort[i] = value;
 392             secondaryErrorShort[i] = error;
 393             secondaryCorrPerNChShort[i] = value/(trackmultiplicityShort[i])/5.0;
 394             secondaryErrorPerNChShort[i] = error/(trackmultiplicityShort[i])/5.0;
 395           }
 396           i++;
 397         }
 398        mysecondaryShortfile.close();
 399     }
 400     //Begin reading in no eff corr corrections
 401     TString secondaryInfileName = "Secondaries"+detector+"NoEffCorr.dat";
 402     ifstream mysecondaryfile3 (secondaryInfileName.Data());
 403     value = 0;
 404     error = 0;
 405     i=0;
 406     if (mysecondaryfile3.is_open()){
 407       while ( mysecondaryfile3.good() )
 408         {
 409           getline (mysecondaryfile3,inline);
 410           istringstream tmp(inline);
 411           tmp >> value;
 412           tmp >> error;
 413           if(i<20){
 414             secondaryCorrNoEffCorr[i] = value;
 415             secondaryErrorNoEffCorr[i] = error;
 416             if(trackmultiplicity[i]>0){
 417               secondaryCorrPerNChNoEffCorr[i] = value/(trackmultiplicity[i])/5.0;
 418               secondaryErrorPerNChNoEffCorr[i] = error/(trackmultiplicity[i])/5.0;
 419             }
 420           }
 421           //cout<<"secondarycorr cb "<<i<<" "<<value<<" +/- "<<error<<" "<<secondaryCorrPerNChNoEffCorr[i]<<" +/- "<<secondaryErrorPerNChNoEffCorr[i]<<endl;
 422           i++;
 423         }
 424         mysecondaryfile3.close();
 425     }
 426     TString secondaryShortInfileName = "Secondaries"+detector+"NoEffCorrShort.dat";
 427     ifstream mysecondaryShortfile3 (secondaryShortInfileName.Data());
 428     i=0;
 429     if (mysecondaryShortfile3.is_open()){
 430       while ( mysecondaryShortfile3.good() && i<10 )
 431         {
 432           getline (mysecondaryShortfile3,inline);
 433           istringstream tmp(inline);
 434           tmp >> value;
 435           tmp >> error;
 436           if(i<10){
 437             secondaryCorrShortNoEffCorr[i] = value;
 438             secondaryErrorShortNoEffCorr[i] = error;
 439             secondaryCorrPerNChShortNoEffCorr[i] = value/(trackmultiplicityShort[i])/5.0;
 440             secondaryErrorPerNChShortNoEffCorr[i] = error/(trackmultiplicityShort[i])/5.0;
 441           }
 442           //cout<<"secondarycorr cb "<<i<<" "<<value<<" +/- "<<error<<" "<<secondaryCorrPerNChShortNoEffCorr[i]<<" +/- "<<secondaryErrorPerNChShortNoEffCorr[i]<<endl;
 443           //cout<<"secondarycorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 444           i++;
 445         }
 446        mysecondaryShortfile3.close();
 447     }
 448
 449 }
 450 void ReadInKaonCorrections(){
 451   cout<<"Reading in kaon corrections..."<<endl;
 452   //junk.PHOS.CutNum6.txt
 453   TString cutstring = "7";
 454   if(detector.Contains("P")){ cutstring = "6";}
 455   if(CutSet==1){
 456     cutstring = "8";
 457   }
 458   if(CutSet==2){
 459     cutstring = "9";
 460   }
 461   TString kaonInfileName = "../spectrafits/KaonCut"+cutstring+"EMCal.dat";
 462   if(detector.Contains("P")){
 463     kaonInfileName = "../spectrafits/KaonCut"+cutstring+"PHOS.dat";
 464   }
 465   cout<<"Reading in "<<kaonInfileName<<endl;
 466     ifstream mykaonfile (kaonInfileName.Data());
 467     string inline;
 468     float value = 0;
 469     float error = 0;
 470     int i=0;
 471     if (mykaonfile.is_open()){
 472       while ( mykaonfile.good() )
 473         {
 474           getline (mykaonfile,inline);
 475           istringstream tmp(inline);
 476           tmp >> value;
 477           tmp >> error;
 478           if(i<10){
 479             kaonCorrShort[i] = value;
 480             kaonErrorShort[i] = error;
 481             kaonCorrPerNChShort[i] = value/(trackmultiplicityShort[i]);
 482             kaonErrorPerNChShort[i] = error/(trackmultiplicityShort[i]);
 483             // cout<<"kaoncorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 484           }
 485           i++;
 486         }
 487         mykaonfile.close();
 488     }
 489
 490
 491   TGraphErrors *graphKaonCorrectionShort = GetKaonCorrectionGraphShort();
 492   int shortbin = 0;
 493   for(int i=0;i<19;i++){
 494     //cout<<"Long bin "<<i<<" short bin "<<shortbin;
 495     if(i<2){//we have exact numbers so we don't need to interpolate
 496       kaonCorr[i] = kaonCorrShort[i];
 497       kaonError[i] = kaonErrorShort[i];
 498       shortbin++;
 499     }
 500     else{
 501       kaonCorr[i] = graphKaonCorrectionShort->Eval(trackmultiplicity[i]) * trackmultiplicity[i];
 502       int altbin = shortbin-1;
 503       if(i%2==1){altbin = shortbin+1;}
 504       //cout<<" altbin "<<altbin;
 505       if(kaonErrorPerNChShort[shortbin]>kaonErrorPerNChShort[altbin]) kaonError[i] = kaonErrorPerNChShort[shortbin] * trackmultiplicity[i];
 506       else{kaonError[i] =  kaonErrorPerNChShort[altbin] * trackmultiplicity[i];}
 507       if(i%2==1 && shortbin<10){shortbin++;}
 508     }
 509     //cout<<"kaoncorr cb "<<i<<" "<<kaonCorr[i]<<" +/- "<<kaonError[i]<<endl;
 510     kaonCorrPerNCh[i] = kaonCorr[i]/(trackmultiplicity[i]);
 511     kaonErrorPerNCh[i] = kaonError[i]/(trackmultiplicity[i]);
 512     //cout<<"min et corr cb "<<i<<" "<< kaonCorr[i]<<" +/- "<<kaonError[i]<<endl;
 513     //cout<<endl;
 514   }
 515   delete graphKaonCorrectionShort;
 516
 517
 518     for(int i=0;i<10;i++){
 519       kaonYieldTotErr[i] = TMath::Sqrt(kaonYieldStatErr[i]*kaonYieldStatErr[i]+kaonYieldSysErr[i]*kaonYieldSysErr[i]);
 520       kaonYieldPerNCh[i] = kaonYield[i]/(trackmultiplicityShort[i])/5;
 521       kaonYieldPerNChErr[i] = kaonYieldTotErr[i]/(trackmultiplicityShort[i])/5;
 522       float total = kaonPlusEt[0][i]+kaonMinusEt[0][i];
 523       float err = kaonPlusEt[1][i]+kaonMinusEt[1][i];
 524       kaonEtPerNCh[i] = total/(trackmultiplicityShort[i])/10;
 525       kaonEtPerNChErr[i] = err/(trackmultiplicityShort[i])/10;
 526     }
 527     //corrections for no eff corr
 528   TString kaonInfileName3 = "../spectrafits/KaonCut"+cutstring+"EMCalNoEffCorr.dat";
 529   if(detector.Contains("P")){
 530     kaonInfileName3 = "../spectrafits/KaonCut"+cutstring+"PHOSNoEffCorr.dat";
 531   }
 532   cout<<"Reading in "<<kaonInfileName3<<endl;
 533   ifstream mykaonfile3 (kaonInfileName3.Data());
 534   value = 0;
 535   error = 0;
 536   i=0;
 537   if (mykaonfile3.is_open()){
 538     while ( mykaonfile3.good() )
 539       {
 540         getline (mykaonfile3,inline);
 541         istringstream tmp(inline);
 542         tmp >> value;
 543         tmp >> error;
 544         if(i<10){
 545           kaonCorrShortNoEffCorr[i] = value;
 546             kaonErrorShortNoEffCorr[i] = error;
 547             kaonCorrPerNChShortNoEffCorr[i] = value/(trackmultiplicityShort[i]);
 548             kaonErrorPerNChShortNoEffCorr[i] = error/(trackmultiplicityShort[i]);
 549             //cout<<"kaoncorr cb "<<i<<" "<<value<<" +/- "<<error<<endl;
 550         }
 551         i++;
 552       }
 553     mykaonfile3.close();
 554   }
 555
 556
 557   TGraphErrors *graphKaonCorrectionShort2 = GetKaonCorrectionGraphShortNoEffCorr();
 558   shortbin = 0;
 559   for(int i=0;i<19;i++){
 560     //cout<<"Long bin "<<i<<" short bin "<<shortbin;
 561     if(i<2){//we have exact numbers so we don't need to interpolate
 562       kaonCorrNoEffCorr[i] = kaonCorrShortNoEffCorr[i];
 563       kaonErrorNoEffCorr[i] = kaonErrorShortNoEffCorr[i];
 564       shortbin++;
 565     }
 566     else{
 567       kaonCorrNoEffCorr[i] = graphKaonCorrectionShort2->Eval(trackmultiplicity[i]) * trackmultiplicity[i];
 568       int altbin = shortbin-1;
 569       if(i%2==1){altbin = shortbin+1;}
 570       //cout<<" altbin "<<altbin;
 571       if(kaonErrorPerNChShortNoEffCorr[shortbin]>kaonErrorPerNChShortNoEffCorr[altbin]) kaonErrorNoEffCorr[i] = kaonErrorPerNChShortNoEffCorr[shortbin] * trackmultiplicity[i];
 572       else{kaonErrorNoEffCorr[i] =  kaonErrorPerNChShortNoEffCorr[altbin] * trackmultiplicity[i];}
 573       if(i%2==1 && shortbin<10){shortbin++;}
 574     }
 575     //cout<<"kaoncorr cb "<<i<<" "<<kaonCorrNoEffCorr[i]<<" +/- "<<kaonErrorNoEffCorr[i]<<endl;
 576     kaonCorrPerNChNoEffCorr[i] = kaonCorrNoEffCorr[i]/(trackmultiplicity[i]);
 577     kaonErrorPerNChNoEffCorr[i] = kaonErrorNoEffCorr[i]/(trackmultiplicity[i]);
 578     //cout<<"min et corr cb "<<i<<" "<< kaonCorrNoEffCorr[i]<<" +/- "<<kaonErrorNoEffCorr[i]<<endl;
 579     //cout<<endl;
 580   }
 581   delete graphKaonCorrectionShort2;
 582
 583
 584
 585
 586 }
 587 void CalculateHadronicCorrectionForOneBin(Int_t centbin1, Int_t centbin2, Bool_t isPhos, Bool_t isOver500MeV, Float_t &correction, Float_t &error, Bool_t effCorr){
 588   //cout<<"cb "<<centbin1<<" - "<<centbin2;
 589     fHistMatchedTracksEvspTvsCentEffCorr->GetZaxis()->SetRange(centbin1,centbin2);
 590     fHistMatchedTracksEvspTvsCentEffCorr500MeV->GetZaxis()->SetRange(centbin1,centbin2);
 591     fHistMatchedTracksEvspTvsCent->GetZaxis()->SetRange(centbin1,centbin2);
 592     TH1D *dataEffCorrTmp = NULL;
 593     TH1D *dataEffCorrTmp2 = NULL;
 594     TH1D *dataTmp = NULL;
 595     TH1D *foundTmp = NULL;
 596     TH1D *notfoundTmp = NULL;
 597     dataTmp = (TH1D*)fHistMatchedTracksEvspTvsCent->Project3D("y");
 598     dataTmp->SetName(Form("dataTmp%i",centbin1));
 599     if(isOver500MeV){
 600       dataEffCorrTmp =(TH1D*) fHistMatchedTracksEvspTvsCentEffCorr500MeV->Project3D("y");
 601       dataEffCorrTmp2 =(TH1D*) fHistMatchedTracksEvspTvsCentEffCorr->Project3D("y");
 602       foundTmp = fHistFoundHadronsvsCent500MeV->ProjectionX(Form("Found%iTmp",centbin1),centbin1,centbin2);
 603       notfoundTmp = fHistNotFoundHadronsvsCent500MeV->ProjectionX(Form("NotFound%iTmp",centbin1),centbin1,centbin2);
 604     }
 605     else{
 606       if(effCorr){
 607         dataEffCorrTmp = (TH1D*)fHistMatchedTracksEvspTvsCentEffCorr->Project3D("y");
 608       }
 609       else{
 610         dataEffCorrTmp = (TH1D*)fHistMatchedTracksEvspTvsCent->Project3D("y");
 611       }
 612       dataEffCorrTmp2 = dataEffCorrTmp;
 613       foundTmp = fHistFoundHadronsvsCent->ProjectionX(Form("Found%iTmp",centbin1),centbin1,centbin2);
 614       notfoundTmp = fHistNotFoundHadronsvsCent->ProjectionX(Form("NotFound%iTmp",centbin1),centbin1,centbin2);
 615     }
 616     float nfound = foundTmp->GetMean();//fHistFoundHadronsvsCent->GetBinContent(bin);
 617     float nnotfound = notfoundTmp->GetMean();//fHistNotFoundHadronsvsCent->GetBinContent(bin);
 618     //cout<<" nfound "<<nfound<<" nnotfound "<<nnotfound<<" total "<<nfound+nnotfound;
 619     float scaleLow = 0;
 620     float scaleHigh = 0;
 621     if(centbin1>=refBin){//for peripheral
 622       scaleHigh = 1.01;
 623       scaleLow = 0.97;
 624     }
 625     else{
 626       float scale1 = 1.0;
 627       float scale2 = 1.0;
 628       float refData = ((TH1D*)data[refBin])->GetMean();
 629       float myData = ((TH1D*)dataTmp)->GetMean();
 630       float refDataEffCorr = ((TH1D*)dataEffCorr[refBin])->GetMean();
 631       float myDataEffCorr = ((TH1D*)dataEffCorrTmp2)->GetMean();
 632       cout<<"data Eff corr "<<myDataEffCorr<<" data eff corr ref "<<refDataEffCorr<<endl;
 633       if(TMath::Abs(myData)>1e-5) scale1 = refData/myData;
 634       if(TMath::Abs(myDataEffCorr)>1e-5) scale2 = refDataEffCorr/myDataEffCorr;
 635       if(scale1<scale2){
 636         scaleLow = 0.99*scale1;
 637         scaleHigh = scale2;
 638       }
 639       else{
 640         scaleLow = 0.99*scale2;
 641         scaleHigh = scale1;
 642       }
 643     }
 644     //    if(averageHadEnergy<0){//if this hasn't been filled yet
 645       Float_t low = 100;
 646       Float_t high = -1;
 647
 648       for(int i=refBinHigh;i<=refBin;i++){
 649         Float_t val = ((TH1D*)dataEffCorr[i])->GetMean();
 650         cout<<" i "<<val;
 651         if(val<low) low = val;
 652         if(val>high) high = val;
 653       }
 654       cout<<endl;
 655       averageHadEnergy =0.97* (low+high)/2.0;
 656       averageHadEnergyError = 0.97*(high-low)/2.0;
 657       cout<<" AVERAGE HAD ENERGY "<<averageHadEnergy<<"+/-"<<averageHadEnergyError<<endl;
 658       //}
 659     float myavg = dataEffCorrTmp->GetMean();
 660     float avg = (scaleLow+scaleHigh)/2.0*myavg;
 661     float err = TMath::Abs((scaleLow-scaleHigh))/2.0*myavg;
 662     cout<<"Nominal value "<<avg<<"+/-"<<err<<endl;
 663     //avg = averageHadEnergy;
 664     //err = averageHadEnergyError;
 665     //cout<<" avg "<<avg<<" +/- "<<err;
 666     if(TMath::Abs(avg)<1e-3){
 667       avg = 1e-3;
 668       cerr<<"WARNING:  ERROR NOT CALCULATED CORRECTLY!!"<<endl;//prevents a crash
 669     }
 670     //factor is the fraction to reduce the track-matched ET by to get the true background ET
 671     //corrfac is the factor to multiply by in order to get the fraction of hadrons leaving deposits which come from low pT
 672     float percentEfficiencyError = 0.04;
 673     //    float  factor = 1-0.04;
 674         float  factor = 1-0.03;
 675 //         float corrfac = 1.275-1;
 676 //     float corrfacerr = 0.059 ;
 677         float corrfac = 0.208938;
 678         float corrfacerr = 0.0357719 ;
 679     float eLowAverage = avg;
 680     float eLowAverageErr = err;
 681     if(isPhos){
 682       factor = 1-0.02;
 683       //factor = 1-0.03;
 684       corrfac = 0.183584;
 685       corrfacerr = 0.0393219;
 686 //       corrfac = 1.300-1;
 687 //       corrfacerr = 0.065;
 688     }
 689     if(isOver500MeV){
 690       eLowAverage = 1.0;
 691       eLowAverageErr = 0.05;
 692       if(isPhos){
 693         //fraction ranges from 5% - 26%
 694         corrfac = ((1.0/0.95 + 1.0/0.74)/2.0 - 1)*eLowAverage/avg;
 695         float corrfacerrtmp = (TMath::Abs((1.0/0.95 - 1.0/0.74)/2.0))*eLowAverage/avg;
 696         corrfacerr = corrfac * TMath::Sqrt(TMath::Power(corrfacerrtmp/corrfac,2)+TMath::Power(eLowAverageErr/eLowAverage,2));
 697       }
 698       else{
 699         //fraction ranges from 5% - 17%
 700         corrfac = ((1.0/0.95 + 1.0/0.83)/2.0 - 1)*eLowAverage/avg;
 701         float corrfacerrtmp = (TMath::Abs((1.0/0.95 - 1.0/0.83)/2.0))*eLowAverage/avg;
 702         corrfacerr = corrfac * TMath::Sqrt(TMath::Power(corrfacerrtmp/corrfac,2)+TMath::Power(eLowAverageErr/eLowAverage,2));
 703       }
 704     }
 705     //the energy from low pT is the fraction of tracks that come from low pT tracks times the average energy of these tracks
 706     float eLow = corrfac * (nfound+nnotfound)*eLowAverage;
 707     float eLowErr = TMath::Sqrt(TMath::Power(corrfacerr*(nfound+nnotfound)*eLowAverage,2)+TMath::Power(eLowAverageErr*corrfac* (nfound+nnotfound),2)+TMath::Power(eLow*percentEfficiencyError,2));//error on the hadronic correction
 708
 709     float eNotFound = nnotfound*avg;
 710     float eNotFoundErr = TMath::Sqrt(TMath::Power(err*nnotfound,2)+TMath::Power(percentEfficiencyError*eNotFound,2));//error on the hadronic correction
 711
 712     float y = (corrfac * (nfound+nnotfound) +nnotfound)*avg;
 713     float finalerr = TMath::Sqrt(TMath::Power(corrfacerr*(nfound+nnotfound)*avg,2)+err*err*(TMath::Power(corrfac* (nfound+nnotfound),2)+nnotfound*nnotfound)+TMath::Power(percentEfficiencyError*y,2));//error on the hadronic correction
 714     correction = y;
 715     error = finalerr;
 716     //cout<<"error "<<finalerr/y<<endl;
 717     //error = 0.0;
 718     //cout<<" corr fac "<<corrfac<<" +/- "<<corrfacerr;
 719     //cout<<" eLow "<<eLow<<" +/- "<<eLowErr<<"("<<eLowErr/y<<") eHigh "<<eNotFound<<" +/- "<<eNotFoundErr<<"("<<eNotFoundErr/y<<") total "<<y<<" +/- "<<finalerr<<"("<<finalerr/y <<")";//<<" frac low "<<eLow/y<<endl;
 720     //cout<<endl;
 721     delete dataEffCorrTmp;
 722     delete foundTmp;
 723     delete notfoundTmp;
 724 }
 725 void CalculateHadronCorrections(Bool_t isPhos){
 726   float plotscale = 5.0;
 727   for(int i=0;i<19;i++){
 728     Float_t correction = 0;
 729     Float_t error = 0;//not above 500 GeV, with eff corr
 730     CalculateHadronicCorrectionForOneBin(i+1,i+1,isPhos,kFALSE,correction,error,kTRUE);
 731     hadCorr[i] = correction;//hadCorrEmcal[i];
 732     hadError[i] = error;//hadErrorEmcal[i];
 733     hadronCorrPerNCh[i] = correction/(trackmultiplicity[i])/plotscale;//hadCorrEmcal[i];
 734     hadronErrorPerNCh[i] = error/(trackmultiplicity[i])/plotscale;//hadErrorEmcal[i];
 735     //cout<<"had cor "<<i<<" "<<correction<<" +/- "<<error<< "  "<<  correction/(trackmultiplicity[i])<< " +/- "<<  error/(trackmultiplicity[i]) <<endl;
 736   }
 737
 738   int j=0;
 739   for(int i=0;i<10;i++){
 740     int centbinlow = i+1;
 741     int centbinhigh = i+1;
 742     if(i<2){//These bins are exactly what they should bin in the 20 bin binning
 743       j++;//i=0 j=0; i=1 j=1
 744       centbinlow = j;
 745       centbinhigh = j;
 746     }
 747     else{
 748       centbinlow = j+1;
 749       centbinhigh = j+2;
 750       j+=2;
 751     }
 752     Float_t correction = 0;
 753     Float_t error = 0;//not above 500 GeV, with eff corr, 10 bins
 754     CalculateHadronicCorrectionForOneBin(centbinlow,centbinhigh,isPhos,kFALSE,correction,error,kTRUE);
 755     hadCorrShort[i] = correction;//hadCorrEmcal[i];
 756     hadErrorShort[i] = error;//hadErrorEmcal[i];
 757     hadronCorrPerNChShort[i] = correction/(trackmultiplicityShort[i])/plotscale;//hadCorrEmcal[i];
 758     hadronErrorPerNChShort[i] = error/(trackmultiplicityShort[i])/plotscale;//hadErrorEmcal[i];
 759     //cout<<"had cor "<<i<<" "<<correction<<" +/- "<<error<< "  "<<  correction/(trackmultiplicityShort[i])<< " +/- "<<  error/(trackmultiplicityShort[i]) <<endl;
 760   }
 761   //No eff corr
 762   for(int i=0;i<19;i++){
 763     Float_t correction = 0;
 764     Float_t error = 0;//not above 500 GeV, without eff corr
 765     CalculateHadronicCorrectionForOneBin(i+1,i+1,isPhos,kFALSE,correction,error,kFALSE);
 766     hadCorrNoEffCorr[i] = correction;//hadCorrEmcalNoEffCorr[i];
 767     hadErrorNoEffCorr[i] = error;//hadErrorEmcalNoEffCorr[i];
 768     hadronCorrPerNChNoEffCorr[i] = correction/(trackmultiplicity[i])/plotscale;//hadCorrEmcalNoEffCorr[i];
 769     hadronErrorPerNChNoEffCorr[i] = error/(trackmultiplicity[i])/plotscale;//hadErrorEmcalNoEffCorr[i];
 770     //cout<<"had cor "<<i<<" "<<correction<<" +/- "<<error<< "  "<<  correction/(trackmultiplicityNoEffCorr[i])<< " +/- "<<  error/(trackmultiplicity[i]) <<endl;
 771   }
 772
 773   int j=0;
 774   for(int i=0;i<10;i++){
 775     int centbinlow = i+1;
 776     int centbinhigh = i+1;
 777     if(i<2){//These bins are exactly what they should bin in the 20 bin binning
 778       j++;//i=0 j=0; i=1 j=1
 779       centbinlow = j;
 780       centbinhigh = j;
 781     }
 782     else{
 783       centbinlow = j+1;
 784       centbinhigh = j+2;
 785       j+=2;
 786     }
 787     Float_t correction = 0;
 788     Float_t error = 0;//not above 500 GeV, without eff corr, 10 bins
 789     CalculateHadronicCorrectionForOneBin(centbinlow,centbinhigh,isPhos,kFALSE,correction,error,kFALSE);
 790     hadCorrShortNoEffCorr[i] = correction;//hadCorrEmcalNoEffCorr[i];
 791     hadErrorShortNoEffCorr[i] = error;//hadErrorEmcalNoEffCorr[i];
 792     hadronCorrPerNChShortNoEffCorr[i] = correction/(trackmultiplicityShort[i])/plotscale;//hadCorrEmcalNoEffCorr[i];
 793     hadronErrorPerNChShortNoEffCorr[i] = error/(trackmultiplicityShort[i])/plotscale;//hadErrorEmcalNoEffCorr[i];
 794     //cout<<"had cor "<<i<<" "<<correction<<" +/- "<<error<< "  "<<  correction/(trackmultiplicityShortNoEffCorr[i])<< " +/- "<<  error/(trackmultiplicityShort[i]) <<endl;
 795   }
 796
 797 }
 798
 799 TGraphErrors *GetSignalFractionGraph(){
 800   //TGraphErrors *gr3 = new TGraphErrors(10,npart,secondaryCorrPerNPartShort,npartErrShort,secondaryErrorPerNPartShort);
 801   TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,signalFraction,trackmultiplicityError,signalFractionError);
 802   //TGraphErrors *gr3 = new TGraphErrors(10,npart,minEtCorrShort,npartErrShort,minEtErrorShort);
 803   SetStyles(gr3,29,TColor::kGreen+3);
 804     return gr3;
 805
 806 }
 807
 808 TGraphErrors *GetMinEtCorrectionGraphShort(){
 809   //TGraphErrors *gr3 = new TGraphErrors(10,npart,secondaryCorrPerNPartShort,npartErrShort,secondaryErrorPerNPartShort);
 810   TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,minEtCorrShort,trackmultiplicityShortError,minEtErrorShort);
 811   //TGraphErrors *gr3 = new TGraphErrors(10,npart,minEtCorrShort,npartErrShort,minEtErrorShort);
 812   SetStyles(gr3,29,TColor::kGreen+3);
 813     return gr3;
 814
 815 }
 816 TGraphErrors *GetMinEtCorrectionGraph(){
 817   //TGraphErrors *gr3 = new TGraphErrors(10,npart,secondaryCorrPerNPartShort,npartErrShort,secondaryErrorPerNPartShort);
 818   TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,minEtCorr,trackmultiplicityError,minEtError);
 819   //TGraphErrors *gr3 = new TGraphErrors(10,npart,minEtCorrShort,npartErrShort,minEtErrorShort);
 820   SetStyles(gr3,30,TColor::kGreen+3);
 821     return gr3;
 822
 823 }
 824
 825
 826 TGraphErrors *GetSecondaryCorrectionGraphShort(){
 827   //TGraphErrors *gr3 = new TGraphErrors(10,npart,secondaryCorrPerNPartShort,npartErrShort,secondaryErrorPerNPartShort);
 828   TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,secondaryCorrPerNChShort,trackmultiplicityShortError,secondaryErrorPerNChShort);
 829   SetStyles(gr3,29,TColor::kGreen+3);
 830     return gr3;
 831
 832 }
 833 TGraphErrors *GetSecondaryCorrectionGraph(){
 834   TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,secondaryCorrPerNCh,trackmultiplicityError,secondaryErrorPerNCh);
 835   SetStyles(gr3,30,TColor::kGreen+3);
 836     return gr3;
 837
 838 }
 839
 840 TGraphErrors *GetSecondaryCorrectionGraphNoEffCorr(){
 841   TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,secondaryCorrPerNChNoEffCorr,trackmultiplicityError,secondaryErrorPerNChNoEffCorr);
 842   SetStyles(gr3,30,TColor::kGreen+3);
 843     return gr3;
 844
 845 }
 846 TGraphErrors *GetSecondaryCorrectionGraphShortNoEffCorr(){
 847   //TGraphErrors *gr3 = new TGraphErrors(10,npart,secondaryCorrPerNPartShort,npartErrShort,secondaryErrorPerNPartShort);
 848   TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,secondaryCorrPerNChShortNoEffCorr,trackmultiplicityShortError,secondaryErrorPerNChShortNoEffCorr);
 849   SetStyles(gr3,29,TColor::kGreen+3);
 850     return gr3;
 851
 852 }
 853 TGraphErrors *GetNeutronCorrectionGraph(){
 854     TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,neutronCorrPerNCh,trackmultiplicityError,neutronErrorPerNCh);
 855     SetStyles(gr3,24,TColor::kBlue);
 856     return gr3;
 857
 858 }
 859 TGraphErrors *GetNeutronCorrectionGraphShort(){
 860   //TGraphErrors *gr3 = new TGraphErrors(10,npart,neutronCorrPerNChShort,npartErrShort,neutronErrorPerNChShort);
 861     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,neutronCorrPerNChShort,trackmultiplicityShortError,neutronErrorPerNChShort);
 862     SetStyles(gr3,20,TColor::kBlue);
 863     return gr3;
 864
 865 }
 866 TGraphErrors *GetNeutronCorrectionGraphNoEffCorr(){
 867     TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,neutronCorrPerNChNoEffCorr,trackmultiplicityError,neutronErrorPerNChNoEffCorr);
 868     SetStyles(gr3,24,TColor::kBlue);
 869     return gr3;
 870
 871 }
 872 TGraphErrors *GetNeutronCorrectionGraphShortNoEffCorr(){
 873   //TGraphErrors *gr3 = new TGraphErrors(10,npart,neutronCorrPerNChShort,npartErrShort,neutronErrorPerNChShort);
 874     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,neutronCorrPerNChShortNoEffCorr,trackmultiplicityShortError,neutronErrorPerNChShortNoEffCorr);
 875     SetStyles(gr3,20,TColor::kBlue);
 876     return gr3;
 877
 878 }
 879 TGraphErrors *GetHadronCorrectionGraph(){
 880     TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,hadronCorrPerNCh,trackmultiplicityError,hadronErrorPerNCh);
 881     SetStyles(gr3,25,1);
 882     return gr3;
 883
 884 }
 885 TGraphErrors *GetHadronCorrectionGraphShort(){
 886   //TGraphErrors *gr3 = new TGraphErrors(10,npart,neutronCorrPerNChShort,xpionerr,neutronErrorPerNChShort);
 887     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,hadronCorrPerNChShort,trackmultiplicityShortError,hadronErrorPerNChShort);
 888     SetStyles(gr3,21,1);
 889     return gr3;
 890
 891 }
 892 TGraphErrors *GetHadronCorrectionGraphNoEffCorr(){
 893     TGraphErrors *gr3 = new TGraphErrors(18,trackmultiplicity,hadronCorrPerNChNoEffCorr,trackmultiplicityError,hadronErrorPerNChNoEffCorr);
 894     SetStyles(gr3,25,1);
 895     return gr3;
 896
 897 }
 898 TGraphErrors *GetHadronCorrectionGraphShortNoEffCorr(){
 899   //TGraphErrors *gr3 = new TGraphErrors(10,npart,neutronCorrPerNChShort,xpionerr,neutronErrorPerNChShort);
 900     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,hadronCorrPerNChShortNoEffCorr,trackmultiplicityShortError,hadronErrorPerNChShortNoEffCorr);
 901     SetStyles(gr3,21,1);
 902     return gr3;
 903
 904 }
 905 TGraphErrors *GetKaonCorrectionGraphShort(){
 906   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 907     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonCorrPerNChShort,trackmultiplicityShortError,kaonErrorPerNChShort);
 908     SetStyles(gr3,33,TColor::kRed);
 909     return gr3;
 910
 911 }
 912 TGraphErrors *GetKaonCorrectionGraph(){
 913   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 914     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,kaonCorrPerNCh,trackmultiplicityError,kaonErrorPerNCh);
 915     SetStyles(gr3,27,TColor::kRed);
 916     return gr3;
 917
 918 }
 919 TGraphErrors *GetKaonCorrectionGraphShortNoEffCorr(){
 920   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 921     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonCorrPerNChShortNoEffCorr,trackmultiplicityShortError,kaonErrorPerNChShortNoEffCorr);
 922     SetStyles(gr3,33,TColor::kRed);
 923     return gr3;
 924
 925 }
 926 TGraphErrors *GetKaonCorrectionGraphNoEffCorr(){
 927   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 928     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,kaonCorrPerNChNoEffCorr,trackmultiplicityError,kaonErrorPerNChNoEffCorr);
 929     SetStyles(gr3,27,TColor::kRed);
 930     return gr3;
 931
 932 }
 933 TGraphErrors *GetKaonGraph(){
 934   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 935     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonYieldPerNCh,trackmultiplicityShortError,kaonYieldPerNChErr);
 936     SetStyles(gr3,33,TColor::kBlue);
 937     return gr3;
 938 }
 939 TGraphErrors *GetKaonEtGraph(){
 940   //TGraphErrors *gr3 = new TGraphErrors(10,xpion,kaonCorrPerNChShort,xpionerr,kaonErrorPerNChShort);
 941     TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonEtPerNCh,trackmultiplicityShortError,kaonEtPerNChErr);
 942     SetStyles(gr3,27,TColor::kBlue);
 943     return gr3;
 944 }
 945
 946
 947 //=====================READ IN DATA===================================
 948 Float_t arrayofzeros[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 949 Float_t trackmultiplicity[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 950 Float_t trackmultiplicityShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 951 Float_t clustermultiplicity[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 952 Float_t clustermultiplicityShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 953 Float_t trackmultiplicityError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 954 Float_t trackmultiplicityShortError[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 955 Float_t clustermultiplicityError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 956 Float_t clustermultiplicityShortError[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 957 Float_t matchedtrackmultiplicity[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 958 Float_t matchedtrackmultiplicityPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 959 Float_t matchedtrackmultiplicityPerNCl[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 960 Float_t notmatchedtrackmultiplicity[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 961 Float_t notmatchedtrackmultiplicityPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 962 Float_t notmatchedtrackmultiplicityPerNCl[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 963 Float_t totaltrackmultiplicityPerNCh[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 964 Float_t totaltrackmultiplicityPerNCl[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 965 Float_t partialCorrEtValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 966 Float_t partialCorrEtValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 967 Float_t partialCorrEtError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 968 Float_t partialCorrEtErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 969 Float_t partialCorrEtPerNChValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 970 Float_t partialCorrEtPerNChValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 971 Float_t partialCorrEtPerNChError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 972 Float_t partialCorrEtPerNChErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 973 Float_t partialCorrEtPerNClValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 974 Float_t partialCorrEtPerNClValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 975 Float_t partialCorrEtPerNClError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 976 Float_t partialCorrEtPerNClErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 977 Float_t partialCorrEtPerNPartPairValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 978 Float_t partialCorrEtPerNPartPairValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 979 Float_t partialCorrEtPerNPartPairError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 980 Float_t partialCorrEtPerNPartPairErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 981 Float_t totalCorrectionPerNPartPairValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 982 Float_t totalCorrectionPerNPartPairValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 983 Float_t totalCorrectionPerNPartPairError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 984 Float_t totalCorrectionPerNPartPairErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 985 Float_t totalCorrectionPerNPartPairValuesNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 986 Float_t totalCorrectionPerNPartPairValuesShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 987 Float_t totalCorrectionPerNPartPairErrorNoEffCorr[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 988 Float_t totalCorrectionPerNPartPairErrorShortNoEffCorr[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 989
 990 Float_t rawEtNoEffCorrValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 991 Float_t rawEtNoEffCorrError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 992 Float_t rawEtNoEffCorrPerNChValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 993 Float_t rawEtNoEffCorrPerNChError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 994 Float_t rawEtNoEffCorrPerNClValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 995 Float_t rawEtNoEffCorrPerNClError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 996 Float_t rawEtAllNoEffCorrValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 997 Float_t rawEtAllNoEffCorrError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 998 Float_t rawEtAllNoEffCorrPerNChValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
 999 Float_t rawEtAllNoEffCorrPerNChError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1000 Float_t rawEtAllNoEffCorrPerNClValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1001 Float_t rawEtAllNoEffCorrPerNClError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1002 Float_t rawEtAllValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1003 Float_t rawEtAllError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1004 Float_t rawEtAllPerNChValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1005 Float_t rawEtAllPerNChError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1006 Float_t rawEtAllPerNClValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1007 Float_t rawEtAllPerNClError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1008
1009 Float_t corrEtValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1010 Float_t corrEtValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1011 Float_t corrEtError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1012 Float_t corrEtErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1013 Float_t corrEtPerNPartPairValues[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1014 Float_t corrEtPerNPartPairValuesShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1015 Float_t corrEtPerNPartPairError[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1016 Float_t corrEtPerNPartPairErrorShort[10] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1017 Float_t corrEtValuesFormulaC[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1018 Float_t corrEtErrorFormulaC[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1019 Float_t corrEtPerNPartPairValuesFormulaC[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1020 Float_t corrEtPerNPartPairErrorFormulaC[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1021 Float_t corrEtValuesFormulaB[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1022 Float_t corrEtErrorFormulaB[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1023 Float_t corrEtPerNPartPairValuesFormulaB[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1024 Float_t corrEtPerNPartPairErrorFormulaB[20] = {0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,0.0};
1025
1026 TH2F *fHistNominalRawEt;
1027 TH2F *fHistNominalNonLinLowEt;
1028 TH2F *fHistNominalNonLinHighEt;
1029 TH2F *fHistNominalEffLowEt;
1030 TH2F *fHistNominalEffHighEt;
1031 TH2F *fHistTotRawEt;
1032 TH2F *fHistTotAllRawEt;
1033 TH2F *fHistTotRawEtNoEffCorr;
1034 TH2F *fHistTotAllRawEtNoEffCorr;
1035 TH2F *fHistTotRawEt500MeV;
1036 TH3F  *fHistMatchedTracksEvspTvsCent;
1037 TH3F  *fHistMatchedTracksEvspTvsCentEffCorr;
1038 TH3F  *fHistMatchedTracksEvspTvsCentEffCorr500MeV;
1039 TH2F *fHistFoundHadronsvsCent;
1040 TH2F *fHistNotFoundHadronsvsCent;
1041 TH2F *fHistFoundHadronsvsCent500MeV;
1042 Int_t refBin = 15;//Reference bin for scaling
1043 Int_t refBinHigh = 8;//Reference bin for scaling
1044 TObjArray data(21);
1045 TObjArray dataFits(21);
1046 TObjArray dataEffCorr(21);
1047 TObjArray dataEffCorr500MeV(21);
1048 TObjArray rawEt(21);
1049 TObjArray rawEtShort(11);
1050 TObjArray rawEtNoEffCorr(21);
1051 TObjArray rawEtAll(21);
1052 TObjArray rawEtAllNoEffCorr(21);
1053
1054 void ReadInData(char *filename,TString det){
1055   cout<<"Reading in data..."<<endl;
1056   Bool_t isPhos = kTRUE;
1057   if(det.Contains("Emc")){isPhos=kFALSE;}
1058
1059   TFile *f = TFile::Open(filename, "READ");
1060     if (!f)
1061     {
1062         std::cerr << "Could not open file: " << filename << std::endl;
1063     }
1064
1065     TList *l = dynamic_cast<TList*>(f->Get("out1"));
1066     if (!l)
1067     {
1068         std::cerr << "Could not get object list from: " << filename << std::endl;
1069     }
1070     TString prefix = "fHistNominal";
1071     fHistNominalRawEt =(TH2F*) l->FindObject((prefix+"RawEt"+det+"Rec").Data());
1072     fHistNominalNonLinLowEt = (TH2F*)l->FindObject((prefix+"NonLinLowEt"+det+"Rec").Data());
1073     fHistNominalNonLinHighEt = (TH2F*)l->FindObject((prefix+"NonLinHighEt"+det+"Rec").Data());
1074     fHistNominalEffLowEt = (TH2F*)l->FindObject((prefix+"EffLowEt"+det+"Rec").Data());
1075     fHistNominalEffHighEt = (TH2F*)l->FindObject((prefix+"EffHighEt"+det+"Rec").Data());
1076     fHistTotRawEt =(TH2F*) l->FindObject("fHistTotRawEtEffCorr");
1077     fHistTotAllRawEt = (TH2F*) l->FindObject("fHistTotAllRawEtEffCorr");
1078     fHistTotRawEtNoEffCorr = (TH2F*)l->FindObject("fHistTotRawEt");
1079     fHistTotAllRawEtNoEffCorr = (TH2F*)l->FindObject("fHistTotAllRawEt");
1080     fHistTotRawEt500MeV =(TH2F*) l->FindObject("fHistTotRawEtEffCorr500MeV");
1081
1082     fHistCentVsNchVsNcl =(TH3F*) l->FindObject("fHistCentVsNchVsNclReco");
1083
1084
1085     fHistMatchedTracksEvspTvsCent =(TH3F*) l->FindObject("fHistMatchedTracksEvspTvsCent");
1086     fHistMatchedTracksEvspTvsCentEffCorr =(TH3F*) l->FindObject("fHistMatchedTracksEvspTvsCentEffTMCorr");
1087     fHistMatchedTracksEvspTvsCentEffCorr500MeV = (TH3F*) l->FindObject("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV");
1088     fHistFoundHadronsvsCent = (TH2F*)l->FindObject("fHistFoundHadronsvsCent");
1089     fHistNotFoundHadronsvsCent = (TH2F*)l->FindObject("fHistNotFoundHadronsvsCent");
1090     fHistFoundHadronsvsCent500MeV = (TH2F*)l->FindObject("fHistFoundHadronsvsCent500MeV");
1091     fHistNotFoundHadronsvsCent500MeV = (TH2F*)l->FindObject("fHistNotFoundHadronsvsCent500MeV");
1092     int nbins = 20;
1093     for(int bin = 1; bin<=nbins;bin++){
1094       fHistMatchedTracksEvspTvsCent->GetZaxis()->SetRange(bin,bin);
1095       data[bin] = fHistMatchedTracksEvspTvsCent->Project3D("y");
1096       fHistMatchedTracksEvspTvsCentEffCorr->GetZaxis()->SetRange(bin,bin);
1097       dataEffCorr[bin] = fHistMatchedTracksEvspTvsCentEffCorr->Project3D("y");
1098       fHistMatchedTracksEvspTvsCentEffCorr500MeV->GetZaxis()->SetRange(bin,bin);
1099       dataEffCorr500MeV[bin] = fHistMatchedTracksEvspTvsCentEffCorr500MeV->Project3D("y");
1100       ((TH1D*)data[bin])->SetName(Form("DataEff%i",bin));
1101       ((TH1D*)dataEffCorr[bin])->SetName(Form("DataEffCorr%i",bin));
1102       ((TH1D*)dataEffCorr500MeV[bin])->SetName(Form("DataEffCorr500MeV%i",bin));
1103
1104       rawEt[bin]= fHistTotRawEt->ProjectionX(Form("RawEt%i",bin),bin,bin);
1105       ((TH1D*)rawEt[bin])->SetName(Form("rawEt%i",bin));
1106       partialCorrEtValues[bin-1] = (Float_t)((TH1D*)rawEt[bin])->GetMean();
1107       partialCorrEtError[bin-1] = (Float_t) ((TH1D*)rawEt[bin])->GetMeanError();
1108       partialCorrEtPerNPartPairValues[bin-1] = partialCorrEtValues[bin-1]/(npart[bin-1])/2.0*10;
1109       partialCorrEtPerNPartPairError[bin-1]  =  partialCorrEtError[bin-1]/(npart[bin-1])/2.0*10;
1110
1111       rawEtNoEffCorr[bin]= fHistTotRawEtNoEffCorr->ProjectionX(Form("RawEtNoEffCorr%i",bin),bin,bin);
1112       ((TH1D*)rawEtNoEffCorr[bin])->SetName(Form("rawEtNoEffCorr%i",bin));
1113       rawEtNoEffCorrValues[bin-1] = (Float_t)((TH1D*)rawEtNoEffCorr[bin])->GetMean();
1114       rawEtNoEffCorrError[bin-1] = (Float_t) ((TH1D*)rawEtNoEffCorr[bin])->GetMeanError();
1115
1116       rawEtAllNoEffCorr[bin]= fHistTotAllRawEtNoEffCorr->ProjectionX(Form("RawEtAllNoEffCorr%i",bin),bin,bin);
1117       ((TH1D*)rawEtAllNoEffCorr[bin])->SetName(Form("rawEtAllNoEffCorr%i",bin));
1118       rawEtAllNoEffCorrValues[bin-1] = (Float_t)((TH1D*)rawEtAllNoEffCorr[bin])->GetMean();
1119       rawEtAllNoEffCorrError[bin-1] = (Float_t) ((TH1D*)rawEtAllNoEffCorr[bin])->GetMeanError();
1120
1121       rawEtAll[bin]= fHistTotAllRawEt->ProjectionX(Form("RawEtAll%i",bin),bin,bin);
1122       ((TH1D*)rawEtAll[bin])->SetName(Form("rawEtAll%i",bin));
1123       rawEtAllValues[bin-1] = (Float_t)((TH1D*)rawEtAll[bin])->GetMean();
1124       rawEtAllError[bin-1] = (Float_t) ((TH1D*)rawEtAll[bin])->GetMeanError();
1125       //cout<<"bin "<<bin<<" "<<partialCorrEtValues[bin-1]<<" "<<rawEtNoEffCorrValues[bin-1]<<" "<<rawEtAllNoEffCorrValues[bin-1]<<" "<<rawEtAllValues[bin-1]<<endl;
1126       //cout<<"bin "<<bin<<" eff corr "<<partialCorrEtValues[bin-1]<<" no eff corr "<<rawEtNoEffCorrValues[bin-1]<<" eff "<< rawEtNoEffCorrValues[bin-1]/partialCorrEtValues[bin-1] <<endl;
1127       if(partialCorrEtValues[bin-1]>0) averageEfficiency[bin-1] = rawEtNoEffCorrValues[bin-1]/partialCorrEtValues[bin-1];
1128
1129       TH1D *temp = fHistNominalRawEt->ProjectionX("temp",bin,bin);
1130       float nominal = temp->GetMean();
1131       //cout<<" Mean "<<temp->GetMean()<<" nbins "<<temp->GetNbinsX()<<endl;
1132       delete temp;
1133       temp = fHistNominalNonLinLowEt->ProjectionX("temp",bin,bin);
1134       float nonlinlow = temp->GetMean();
1135       delete temp;
1136       temp = fHistNominalNonLinHighEt->ProjectionX("temp",bin,bin);
1137       float nonlinhigh = temp->GetMean();
1138       delete temp;
1139       temp = fHistNominalEffLowEt->ProjectionX("temp",bin,bin);
1140       float efflow = temp->GetMean();
1141       delete temp;
1142       temp = fHistNominalEffHighEt->ProjectionX("temp",bin,bin);
1143       float effhigh = temp->GetMean();
1144       delete temp;
1145       float nonlinfracerr = 0;
1146       if(nonlinhigh >0 || nonlinlow >0) nonlinfracerr = TMath::Abs(nonlinhigh-nonlinlow)/(nonlinhigh+nonlinlow);
1147       float efffracerr = 0;
1148       if(effhigh >0 || efflow>0)efffracerr = TMath::Abs(effhigh-efflow)/(effhigh+efflow);
1149       //cout<<"cb "<<bin-1<<" nonlinerr "<<nonlinfracerr<<" efficiencyError "<<efffracerr<<endl;
1150       nonLinError[bin-1] = nonlinfracerr;
1151       efficiencyError[bin-1] = efffracerr;
1152
1153       //(not)matchedtrackmultiplicity, (not)matchedtrackmultiplicityPerNCh, (not)matchedtrackmultiplicityPerNCl
1154       temp = fHistFoundHadronsvsCent->ProjectionX(Form("Found%iTmp",bin),bin,bin);
1155       matchedtrackmultiplicity[bin-1] = temp->GetMean();
1156       delete temp;
1157       temp = fHistNotFoundHadronsvsCent->ProjectionX(Form("NotFound%iTmp",bin),bin,bin);
1158       notmatchedtrackmultiplicity[bin-1] = temp->GetMean();
1159       delete temp;
1160
1161
1162
1163     }
1164
1165
1166     for(int cb=0;cb<20;cb++){
1167       fHistCentVsNchVsNcl->GetXaxis()->SetRange(cb+1,cb+1);
1168       TH1D *trackmultiplicityHist = fHistCentVsNchVsNcl->Project3D("y");
1169       TH1D *clustermultiplicityHist = fHistCentVsNchVsNcl->Project3D("z");
1170       trackmultiplicity[cb] = (Float_t) trackmultiplicityHist->GetMean();
1171       clustermultiplicity[cb] = (Float_t)clustermultiplicityHist->GetMean();
1172       trackmultiplicityError[cb] = (Float_t) trackmultiplicityHist->GetMeanError();
1173       clustermultiplicityError[cb] = (Float_t)clustermultiplicityHist->GetMeanError();
1174       delete trackmultiplicityHist;
1175       delete clustermultiplicityHist;
1176     }
1177     int cb1 = 0;
1178     for(int cb=0;cb<10;cb++){
1179       int cb2 = cb1+1;
1180       if(cb1<2) cb2 = cb1;
1181       //cout<<"From "<<cb1<<" to "<<cb2<<endl;
1182       fHistCentVsNchVsNcl->GetXaxis()->SetRange(cb1+1,cb2+1);
1183       TH1D *trackmultiplicityHistShort = fHistCentVsNchVsNcl->Project3D("y");
1184       TH1D *clustermultiplicityHistShort = fHistCentVsNchVsNcl->Project3D("z");
1185       trackmultiplicityShort[cb] = (Float_t) trackmultiplicityHistShort->GetMean();
1186       clustermultiplicityShort[cb] = (Float_t)clustermultiplicityHistShort->GetMean();
1187       trackmultiplicityShortError[cb] = (Float_t) trackmultiplicityHistShort->GetMeanError();
1188       clustermultiplicityShortError[cb] = (Float_t)clustermultiplicityHistShort->GetMeanError();
1189       delete trackmultiplicityHistShort;
1190       delete clustermultiplicityHistShort;
1191       if(cb1<2) cb1++;
1192       else{cb1+=2;}
1193       rawEtShort[cb]= fHistTotRawEt->ProjectionX(Form("RawEtShort%i",bin),cb1+1,cb2+1);
1194       ((TH1D*)rawEtShort[cb])->SetName(Form("rawEtShort%i",cb));
1195       partialCorrEtValuesShort[cb] = (Float_t)((TH1D*)rawEtShort[cb])->GetMean();
1196       partialCorrEtErrorShort[cb] = (Float_t) ((TH1D*)rawEtShort[cb])->GetMeanError();
1197
1198       partialCorrEtPerNChValuesShort[cb] = partialCorrEtValuesShort[cb]/(trackmultiplicityShort[cb])/2.0;
1199       partialCorrEtPerNChErrorShort[cb]  =  partialCorrEtErrorShort[cb]/(trackmultiplicityShort[cb])/2.0;
1200
1201       partialCorrEtPerNPartPairValuesShort[cb] = partialCorrEtValues[cb]/(npart[cb])/2.0*10;
1202       partialCorrEtPerNPartPairErrorShort[cb]  =  partialCorrEtError[cb]/(npart[cb])/2.0*10;
1203
1204       TH1D *temp = fHistNominalRawEt->ProjectionX("temp",cb1+1,cb2+1);
1205       float nominal = temp->GetMean();
1206       //cout<<" Mean "<<temp->GetMean()<<" nbins "<<temp->GetNbinsX()<<endl;
1207       delete temp;
1208       temp = fHistNominalNonLinLowEt->ProjectionX("temp",cb1+1,cb2+1);
1209       float nonlinlow = temp->GetMean();
1210       delete temp;
1211       temp = fHistNominalNonLinHighEt->ProjectionX("temp",cb1+1,cb2+1);
1212       float nonlinhigh = temp->GetMean();
1213       delete temp;
1214       temp = fHistNominalEffLowEt->ProjectionX("temp",cb1+1,cb2+1);
1215       float efflow = temp->GetMean();
1216       delete temp;
1217       temp = fHistNominalEffHighEt->ProjectionX("temp",cb1+1,cb2+1);
1218       float effhigh = temp->GetMean();
1219       delete temp;
1220       float nonlinfracerr = 0;
1221       if(nonlinhigh >0 || nonlinlow >0) nonlinfracerr = TMath::Abs(nonlinhigh-nonlinlow)/(nonlinhigh+nonlinlow);
1222       if(isPhos){
1223         nonlinfracerr = 0.005;
1224       }
1225       float efffracerr = 0;
1226       if(effhigh >0 || efflow>0)efffracerr = TMath::Abs(effhigh-efflow)/(effhigh+efflow);
1227       nonLinErrorShort[cb] = nonlinfracerr;
1228       if(isPhos){
1229         efficiencyErrorShort[cb] = 0.005;
1230       }
1231       else{
1232         efficiencyErrorShort[cb] = 0.02;
1233       }
1234
1235
1236
1237     }
1238
1239 }
1240 void ApplyCorrections(Float_t scale){//scale takes into account the acceptance in eta and phi and the 1/etaacc
1241   //correlation of errors:
1242     //hadCorr - calculated using matched tracks, correcting with tracking efficiency, systematic error dominated by uncertainty in energy deposited in calorimeter
1243   //kaon correction - calculated using kaon spectra, systematic errors from spectra systematic errors
1244   //neutron correction - systematic errors a bit of a fudge.
1245   //secondary correction - systematic errors from Nch vs Ncl scaling
1246   //etmin - calculated from kinematics, simulation, and pion spectra
1247   //efficiency error - determined from different material budgets
1248   //nonlinearity error - difference between test beam data and simulation
1249   //arguably kaon correction and etmin correction are somewhat correlated
1250
1251   for(int cb = 0;cb<19;cb++){
1252     if(trackmultiplicity[cb]>0){
1253       partialCorrEtPerNChValues[cb] = partialCorrEtValues[cb]/(trackmultiplicity[cb]);
1254       partialCorrEtPerNChError[cb]  =  partialCorrEtError[cb]/(trackmultiplicity[cb]);
1255       rawEtNoEffCorrPerNChValues[cb] = rawEtNoEffCorrValues[cb]/(trackmultiplicity[cb]);
1256       rawEtNoEffCorrPerNChError[cb]  =  rawEtNoEffCorrError[cb]/(trackmultiplicity[cb]);
1257       rawEtAllNoEffCorrPerNChValues[cb] = rawEtAllNoEffCorrValues[cb]/(trackmultiplicity[cb]);
1258       rawEtAllNoEffCorrPerNChError[cb]  =  rawEtAllNoEffCorrError[cb]/(trackmultiplicity[cb]);
1259       rawEtAllPerNChValues[cb] = rawEtAllValues[cb]/(trackmultiplicity[cb]);
1260       rawEtAllPerNChError[cb]  =  rawEtAllError[cb]/(trackmultiplicity[cb]);
1261       matchedtrackmultiplicityPerNCh[cb] = matchedtrackmultiplicity[cb]/(trackmultiplicity[cb]);
1262       notmatchedtrackmultiplicityPerNCh[cb] = notmatchedtrackmultiplicity[cb]/(trackmultiplicity[cb]);
1263     }
1264     if(clustermultiplicity[cb]>0){
1265       partialCorrEtPerNClValues[cb] = partialCorrEtValues[cb]/(clustermultiplicity[cb]);
1266       partialCorrEtPerNClError[cb]  =  partialCorrEtError[cb]/(clustermultiplicity[cb]);
1267       rawEtNoEffCorrPerNClValues[cb] = rawEtNoEffCorrValues[cb]/(clustermultiplicity[cb]);
1268       rawEtNoEffCorrPerNClError[cb]  =  rawEtNoEffCorrError[cb]/(clustermultiplicity[cb]);
1269       rawEtAllNoEffCorrPerNClValues[cb] = rawEtAllNoEffCorrValues[cb]/(clustermultiplicity[cb]);
1270       rawEtAllNoEffCorrPerNClError[cb]  =  rawEtAllNoEffCorrError[cb]/(clustermultiplicity[cb]);
1271       rawEtAllPerNClValues[cb] = rawEtAllValues[cb]/(clustermultiplicity[cb]);
1272       rawEtAllPerNClError[cb]  =  rawEtAllError[cb]/(clustermultiplicity[cb]);
1273       matchedtrackmultiplicityPerNCl[cb] = matchedtrackmultiplicity[cb]/(clustermultiplicity[cb]);
1274       notmatchedtrackmultiplicityPerNCl[cb] = notmatchedtrackmultiplicity[cb]/(clustermultiplicity[cb]);
1275     }
1276
1277     totaltrackmultiplicityPerNCh[cb] = matchedtrackmultiplicityPerNCh[cb] + notmatchedtrackmultiplicityPerNCh[cb];
1278     totaltrackmultiplicityPerNCl[cb] = matchedtrackmultiplicityPerNCl[cb] + notmatchedtrackmultiplicityPerNCl[cb];
1279
1280     //cout<<"cb "<<cb<<" "<<partialCorrEtValues[cb]<<" "<<rawEtNoEffCorrValues[cb]<<" "<<rawEtAllNoEffCorrValues[cb]<<" "<<rawEtAllValues[cb]<<endl;
1281     //cout<<"cb "<<cb<<" "<<partialCorrEtPerNChValues[cb]<<" "<<rawEtNoEffCorrPerNChValues[cb]<<" "<<rawEtAllNoEffCorrPerNChValues[cb]<<" "<<rawEtAllPerNChValues[cb]<<endl;
1282       //cout<<"cb "<<cb<<" "<<partialCorrEtPerNChValues[cb]<<" "<< partialCorrEtValues[cb]<<" "<< partialCorrEtError[cb]<<" "<<trackmultiplicity[cb] <<endl;
1283     corrEtValues[cb] = scale*(partialCorrEtValues[cb] - hadCorr[cb] - kaonCorr[cb] - neutronCorr[cb] - secondaryCorr[cb])/minEtCorr[cb];
1284     //cout<<"cb "<<cb<<"\t"<<corrEtValues[cb]<<" = \t"<<scale<<"*\t("<<partialCorrEtValues[cb]<<" -\t"<<hadCorr[cb]<<" -\t"<<kaonCorr[cb]<<" -\t"<<neutronCorr[cb]<<" -\t"<<secondaryCorr[cb]<<")\t/"<<minEtCorr[cb]<<endl;
1285     corrEtValuesFormulaB[cb] = scale*(rawEtNoEffCorrValues[cb] - hadCorrNoEffCorr[cb] - kaonCorrNoEffCorr[cb] - neutronCorrNoEffCorr[cb] - secondaryCorrNoEffCorr[cb])/minEtCorr[cb]/averageEfficiency[cb];
1286     //cout<<"cb "<<corrEtValuesFormulaB[cb]<<" = "<<scale<<"*("<<rawEtNoEffCorrValues[cb]<<" - "<<hadCorrNoEffCorr[cb]<<" - "<<kaonCorrNoEffCorr[cb]<<" - "<<neutronCorrNoEffCorr[cb]<<" - "<<secondaryCorrNoEffCorr[cb]<<")/"<<minEtCorr[cb]<<"/"<<averageEfficiency[cb]<<endl;
1287     //cout<<" fractions: had "<< hadCorrNoEffCorr[cb]/rawEtNoEffCorrValues[cb] <<" kaon "<<kaonCorrNoEffCorr[cb]/rawEtNoEffCorrValues[cb]<<" neutron "<<neutronCorrNoEffCorr[cb]/rawEtNoEffCorrValues[cb]<<" secondary "<<secondaryCorrNoEffCorr[cb]/rawEtNoEffCorrValues[cb]<<endl;
1288     //signalFraction[cb] = (1.0 - (hadCorrNoEffCorr[cb] + kaonCorrNoEffCorr[cb] + neutronCorrNoEffCorr[cb] + secondaryCorrNoEffCorr[cb])/rawEtNoEffCorrValues[cb]);
1289     signalFraction[cb] = (1.0 - (hadCorr[cb] + kaonCorr[cb] + neutronCorr[cb] + secondaryCorr[cb])/partialCorrEtValues[cb]);
1290     //signalFraction[cb] = 0.25;
1291     //cout<<"cb "<<cb<<" fractions: \thad "<< hadCorr[cb]/partialCorrEtValues[cb]<<"+/-" << hadError[cb]/partialCorrEtValues[cb]<<"\tkaon "<<kaonCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<kaonError[cb]/partialCorrEtValues[cb]<<"\tneutron "<<neutronCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<neutronError[cb]/partialCorrEtValues[cb]<<"\tsecondary "<<secondaryCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<secondaryError[cb]/partialCorrEtValues[cb]<<endl;//rawEtValues
1292     corrEtValuesFormulaC[cb] = scale*partialCorrEtValues[cb]*signalFraction[cb]/minEtCorr[cb];
1293     //corrEtValuesFormulaC[cb] = scale*rawEtNoEffCorrValues[cb]*signalFraction[cb]/minEtCorr[cb]/0.22;
1294     //cout<<"cb "<<cb<<" "<<corrEtValuesFormulaC[cb]<<" = "<<scale<<"*"<<rawEtNoEffCorrValues[cb]<<"*"<<signalFraction[cb]<<"/"<<minEtCorr[cb]<<"/"<<averageEfficiency[cb]<<endl;
1295
1296     totalCorrectionPerNPartPairValues[cb] = hadCorr[cb] + kaonCorr[cb] + neutronCorr[cb] + secondaryCorr[cb];
1297     totalCorrectionPerNPartPairValuesNoEffCorr[cb] = hadCorrNoEffCorr[cb] + kaonCorrNoEffCorr[cb] + neutronCorrNoEffCorr[cb] + secondaryCorrNoEffCorr[cb];
1298     //this comes up enough in the error calculations we'll just make a variable for it
1299     float  partialEt = scale*(partialCorrEtValues[cb])/minEtCorr[cb];
1300     float  partialEtNoEffCorr = scale*(rawEtNoEffCorrValues[cb])/minEtCorr[cb];
1301     //add up the error squared
1302     float err = 0;
1303     float partialerr = 0;
1304     float totalcorrpartialerr = 0;
1305     float errNoEffCorr = 0;
1306     float partialerrNoEffCorr = 0;
1307     float totalcorrpartialerrNoEffCorr = 0;
1308     bool writeerror = false;
1309     if(writeerror){
1310       cout<<"fraction errors: min et "<<minEtError[cb]/minEtCorr[cb];
1311     }
1312     //if(writeerror)cout<<"partialEt "<<partialEt<<" err^2 = ";
1313     //Et min correction
1314     //partialerr += TMath::Power(minEtError[cb]/minEtCorr[cb]*corrEtValues[cb],2);
1315     partialerr = TMath::Power(minEtError[cb]/minEtCorr[cb]*corrEtValues[cb],2.0);
1316     partialerrNoEffCorr = TMath::Power(minEtError[cb]/minEtCorr[cb]*corrEtValuesFormulaB[cb],2.0);
1317     //if(writeerror)cout<<partialerr<<"+";
1318     err+=partialerr;
1319     errNoEffCorr+=partialerrNoEffCorr;
1320     //nonlinearity correction - this is saved as a fractional error
1321     partialerr = TMath::Power(nonLinError[cb]*corrEtValues[cb],2.0);
1322     partialerrNoEffCorr = TMath::Power(nonLinError[cb]*corrEtValuesFormulaB[cb],2.0);
1323     if(writeerror){cout<<" nonlinearity "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1324     //if(writeerror)cout<<partialerr<<"+";
1325     err+=partialerr;
1326     errNoEffCorr+=partialerrNoEffCorr;
1327     //efficiency correction - this is also saved as a fractional error
1328     partialerr = TMath::Power(efficiencyError[cb]*corrEtValues[cb],2);
1329     partialerrNoEffCorr = TMath::Power(efficiencyError[cb]*corrEtValues[cb],2);
1330     if(writeerror){cout<<" efficiency "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1331     //if(writeerror)cout<<partialerr<<"+";
1332     err+=partialerr;
1333     errNoEffCorr+=partialerrNoEffCorr;
1334     //hadron correction
1335     partialerr = TMath::Power(hadError[cb]*scale/minEtCorr[cb],2);
1336     partialerrNoEffCorr = TMath::Power(hadErrorNoEffCorr[cb]*scale/minEtCorr[cb],2);
1337     totalcorrpartialerr += TMath::Power(hadError[cb],2);
1338     totalcorrpartialerrNoEffCorr += TMath::Power(hadErrorNoEffCorr[cb],2);
1339     if(writeerror){cout<<" hadronic corr "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1340     //if(writeerror)cout<<partialerr<<"+";
1341     err+=partialerr;
1342     errNoEffCorr+=partialerrNoEffCorr;
1343     //neutron correction
1344     partialerr = TMath::Power(neutronError[cb]*scale/minEtCorr[cb],2);
1345     totalcorrpartialerr += TMath::Power(neutronError[cb],2);
1346     partialerrNoEffCorr = TMath::Power(neutronErrorNoEffCorr[cb]*scale/minEtCorr[cb],2);
1347     totalcorrpartialerrNoEffCorr += TMath::Power(neutronErrorNoEffCorr[cb],2);
1348     if(writeerror){cout<<" neutron corr "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1349     //if(writeerror)cout<<partialerr<<"+";
1350     err+=partialerr;
1351     errNoEffCorr+=partialerrNoEffCorr;
1352     //kaon correction
1353     partialerr = TMath::Power(kaonError[cb]*scale/minEtCorr[cb],2);
1354     totalcorrpartialerr += TMath::Power(kaonError[cb],2);
1355     partialerrNoEffCorr = TMath::Power(kaonErrorNoEffCorr[cb]*scale/minEtCorr[cb],2);
1356     totalcorrpartialerrNoEffCorr += TMath::Power(kaonErrorNoEffCorr[cb],2);
1357     if(writeerror){cout<<" kaon corr "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1358     //if(writeerror)cout<<partialerr<<"+";
1359     err+=partialerr;
1360     errNoEffCorr+=partialerrNoEffCorr;
1361     //secondary correction
1362     partialerr = TMath::Power(secondaryError[cb]*scale/minEtCorr[cb],2);
1363     totalcorrpartialerr += TMath::Power(secondaryError[cb],2);
1364     partialerrNoEffCorr = TMath::Power(secondaryErrorNoEffCorr[cb]*scale/minEtCorr[cb],2);
1365     totalcorrpartialerrNoEffCorr += TMath::Power(secondaryErrorNoEffCorr[cb],2);
1366     if(writeerror){cout<<" secondaries "<<TMath::Sqrt(partialerr)/corrEtValues[cb];}
1367     //if(writeerror)cout<<partialerr;
1368     err+=partialerr;
1369     errNoEffCorr+=partialerrNoEffCorr;
1370     //And take the square root
1371     err = TMath::Sqrt(err);
1372     errNoEffCorr = TMath::Sqrt(errNoEffCorr);
1373     totalCorrectionPerNPartPairError[cb] = TMath::Sqrt(totalcorrpartialerr);
1374     totalCorrectionPerNPartPairErrorNoEffCorr[cb] = TMath::Sqrt(totalcorrpartialerrNoEffCorr);
1375     if(writeerror)cout<<" = "<<err/corrEtValues[cb]<<endl;
1376     signalFractionError[cb] = totalCorrectionPerNPartPairError[cb]/partialCorrEtValues[cb];
1377     //cout<<"signal fraction "<<signalFraction[cb]<<" +/- ";
1378     //cout<<"cb "<<cb<<" fractions: \thad "<< hadCorr[cb]/partialCorrEtValues[cb]<<"+/-" << hadError[cb]/partialCorrEtValues[cb]<<"\tkaon "<<kaonCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<kaonError[cb]/partialCorrEtValues[cb]<<"\tneutron "<<neutronCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<neutronError[cb]/partialCorrEtValues[cb]<<"\tsecondary "<<secondaryCorr[cb]/partialCorrEtValues[cb]<<"+/-"<<secondaryError[cb]/partialCorrEtValues[cb]<<"\tsignalfrac "<<signalFraction[cb]<<"+/-"<<signalFractionError[cb]<<endl;//rawEtValues
1379     if(partialCorrEtValues[cb]>0)cout<<totalCorrectionPerNPartPairError[cb]/partialCorrEtValues[cb];
1380     //cout<<endl;
1381     corrEtError[cb] = err;
1382     corrEtErrorFormulaB[cb] = errNoEffCorr;
1383     corrEtPerNPartPairValues[cb] = corrEtValues[cb]/(npart[cb]/2.0);
1384     corrEtPerNPartPairError[cb]  =  corrEtError[cb]/(npart[cb]/2.0);
1385     corrEtPerNPartPairValuesFormulaC[cb] = corrEtValuesFormulaC[cb]/(npart[cb]/2.0);
1386     corrEtPerNPartPairErrorFormulaC[cb]  =  corrEtErrorFormulaC[cb]/(npart[cb]/2.0);
1387     corrEtPerNPartPairValuesFormulaB[cb] = corrEtValuesFormulaB[cb]/(npart[cb]/2.0);
1388     corrEtPerNPartPairErrorFormulaB[cb]  =  corrEtErrorFormulaB[cb]/(npart[cb]/2.0);
1389     totalCorrectionPerNPartPairValues[cb] = totalCorrectionPerNPartPairValues[cb]/20.0/trackmultiplicity[cb];
1390     totalCorrectionPerNPartPairError[cb] = totalCorrectionPerNPartPairError[cb]/20.0/trackmultiplicity[cb];
1391     totalCorrectionPerNPartPairValuesNoEffCorr[cb] = totalCorrectionPerNPartPairValuesNoEffCorr[cb]/20.0/trackmultiplicity[cb];
1392     totalCorrectionPerNPartPairErrorNoEffCorr[cb] = totalCorrectionPerNPartPairErrorNoEffCorr[cb]/20.0/trackmultiplicity[cb];
1393     //cout<<"test cb "<<cb<<" total corr/npart pair "<<totalCorrectionPerNPartPairValues[cb]<<" "<<totalCorrectionPerNPartPairError[cb]<<endl;
1394     //cout<<"cb "<<cb <<" et "<< corrEtPerNPartPairValues[cb] <<" +/- "<<corrEtPerNPartPairError[cb];
1395     //cout<<"cb "<<cb <<" et "<< corrEtPerNPartPairValuesFormulaC[cb] <<" +/- "<<corrEtPerNPartPairErrorFormulaC[cb]<<endl;
1396
1397     //cout<<" = "<<scale<<"*"<<"("<<partialCorrEtValues[cb]<<"+/-"<<partialCorrEtError[cb]<<" - "<<hadCorr[cb]<<"+/-"<<hadError[cb]<<" - "<<kaonCorr[cb]<<"+/-"<<kaonError[cb]<<" - "<<neutronCorr[cb]<<"+/-"<<neutronError[cb]<<" - "<<secondaryCorr[cb]<<"+/-"<<secondaryError[cb]<<")/"<<minEtCorr[cb]<<"+/-"<<minEtError[cb];
1398     //cout<<endl;
1399   }
1400 //   for(int cb = 0;cb<10;cb++){
1401 //   }
1402
1403 }
1404
1405 TGraphErrors *GetEtGraph(){
1406     TGraphErrors *gr3 = new TGraphErrors(20,npart,corrEtPerNPartPairValues,npartErr,corrEtPerNPartPairError);
1407     SetStyles(gr3,25,1);
1408     return gr3;
1409 }
1410
1411 TGraphErrors *GetEtGraphFormulaC(){
1412     TGraphErrors *gr3 = new TGraphErrors(20,npart,corrEtPerNPartPairValuesFormulaC,npartErr,corrEtPerNPartPairErrorFormulaC);
1413     for(int i=0;i<20;i++){
1414       //cout<<"i "<<i<<" "<<clustermultiplicity[i]<<": "<<rawEtNoEffCorrPerNClValues[i]<<"+/-"<<rawEtNoEffCorrPerNClError[i]<<endl;
1415     }
1416     SetStyles(gr3,21,1);
1417     return gr3;
1418 }
1419 TGraphErrors *GetEtGraphFormulaB(){
1420     TGraphErrors *gr3 = new TGraphErrors(20,npart,corrEtPerNPartPairValuesFormulaB,npartErr,corrEtPerNPartPairErrorFormulaB);
1421     for(int i=0;i<20;i++){
1422       //cout<<"i "<<i<<" "<<clustermultiplicity[i]<<": "<<rawEtNoEffCorrPerNClValues[i]<<"+/-"<<rawEtNoEffCorrPerNClError[i]<<endl;
1423     }
1424     SetStyles(gr3,29,1);
1425     return gr3;
1426 }
1427 TGraphErrors *GetPartialCorrEtPerNChGraph(){
1428     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,partialCorrEtPerNChValues,trackmultiplicityError,partialCorrEtPerNChError);
1429     //TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonEtPerNCh,trackmultiplicityShortError,kaonEtPerNChErr);
1430
1431     SetStyles(gr3,25,1);
1432     return gr3;
1433 }
1434 TGraphErrors *GetPartialCorrEtPerNClGraph(){
1435     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,partialCorrEtPerNClValues,clustermultiplicityError,partialCorrEtPerNClError);
1436     //TGraphErrors *gr3 = new TGraphErrors(10,trackmultiplicityShort,kaonEtPerNCh,trackmultiplicityShortError,kaonEtPerNChErr);
1437
1438     SetStyles(gr3,25,1);
1439     return gr3;
1440 }
1441 TGraphErrors *GetRawEtNoEffCorrPerNChGraph(){
1442     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,rawEtNoEffCorrPerNChValues,trackmultiplicityError,rawEtNoEffCorrPerNChError);
1443     SetStyles(gr3,21,1);
1444     return gr3;
1445 }
1446 TGraphErrors *GetRawEtNoEffCorrPerNClGraph(){
1447     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,rawEtNoEffCorrPerNClValues,clustermultiplicityError,rawEtNoEffCorrPerNClError);
1448     for(int i=0;i<20;i++){
1449       //cout<<"i "<<i<<" "<<clustermultiplicity[i]<<": "<<rawEtNoEffCorrPerNClValues[i]<<"+/-"<<rawEtNoEffCorrPerNClError[i]<<endl;
1450     }
1451     SetStyles(gr3,21,1);
1452     return gr3;
1453 }
1454 TGraphErrors *GetRawEtAllNoEffCorrPerNChGraph(){
1455     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,rawEtAllNoEffCorrPerNChValues,trackmultiplicityError,rawEtAllNoEffCorrPerNChError);
1456     SetStyles(gr3,26,1);
1457     return gr3;
1458 }
1459 TGraphErrors *GetRawEtAllNoEffCorrPerNClGraph(){
1460     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,rawEtAllNoEffCorrPerNClValues,clustermultiplicityError,rawEtAllNoEffCorrPerNClError);
1461     SetStyles(gr3,26,1);
1462     return gr3;
1463 }
1464 TGraphErrors *GetRawEtAllPerNChGraph(){
1465     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,rawEtAllPerNChValues,trackmultiplicityError,rawEtAllPerNChError);
1466     SetStyles(gr3,26,1);
1467     return gr3;
1468 }
1469 TGraphErrors *GetRawEtAllPerNClGraph(){
1470     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,rawEtAllPerNClValues,clustermultiplicityError,rawEtAllPerNClError);
1471     SetStyles(gr3,22,1);
1472     return gr3;
1473 }
1474 TGraphErrors *GetMatchedTracksPerNChGraph(){
1475     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,matchedtrackmultiplicityPerNCh,trackmultiplicityError,arrayofzeros);
1476     SetStyles(gr3,20,1);
1477     return gr3;
1478 }
1479 TGraphErrors *GetMatchedTracksPerNClGraph(){
1480     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,matchedtrackmultiplicityPerNCl,clustermultiplicityError,arrayofzeros);
1481     SetStyles(gr3,20,1);
1482     return gr3;
1483 }
1484 TGraphErrors *GetTotalTracksPerNClGraph(){
1485     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,totaltrackmultiplicityPerNCl,clustermultiplicityError,arrayofzeros);
1486     SetStyles(gr3,34,1);
1487     return gr3;
1488 }
1489 TGraphErrors *GetTotalTracksPerNChGraph(){
1490     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,totaltrackmultiplicityPerNCh,trackmultiplicityError,arrayofzeros);
1491     SetStyles(gr3,34,1);
1492     return gr3;
1493 }
1494 TGraphErrors *GetNotMatchedTracksPerNChGraph(){
1495     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,notmatchedtrackmultiplicityPerNCh,trackmultiplicityError,arrayofzeros);
1496     SetStyles(gr3,24,1);
1497     return gr3;
1498 }
1499 TGraphErrors *GetNotMatchedTracksPerNClGraph(){
1500     TGraphErrors *gr3 = new TGraphErrors(20,clustermultiplicity,notmatchedtrackmultiplicityPerNCl,clustermultiplicityError,arrayofzeros);
1501     SetStyles(gr3,24,1);
1502     return gr3;
1503 }
1504 TGraphErrors *GetPartialCorrEtPerNPartPairGraph(){
1505     TGraphErrors *gr3 = new TGraphErrors(20,npart,partialCorrEtPerNPartPairValues,npartErr,partialCorrEtPerNPartPairError);
1506     SetStyles(gr3,22,1);
1507     return gr3;
1508 }
1509 TGraphErrors *GetTotalCorrectionPerNChGraph(){
1510     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,totalCorrectionPerNPartPairValues,trackmultiplicityError,totalCorrectionPerNPartPairError);
1511     SetStyles(gr3,34,TColor::kOrange);
1512     return gr3;
1513 }
1514 TGraphErrors *GetTotalCorrectionPerNChGraphNoEffCorr(){
1515     TGraphErrors *gr3 = new TGraphErrors(20,trackmultiplicity,totalCorrectionPerNPartPairValuesNoEffCorr,trackmultiplicityError,totalCorrectionPerNPartPairErrorNoEffCorr);
1516     SetStyles(gr3,34,TColor::kOrange);
1517     return gr3;
1518 }
1519 //=================================Plotting code====================================
1520 Bool_t sim = false;
1521 //Bool_t isPhos = kFALSE;
1522 TString detector = "";
1523 //void PlotEmEtVer2(TString filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.root")
1524 void PlotEmEtVer2(Bool_t isPhos = kFALSE, Bool_t isMC = kFALSE, Int_t cutset = 0)
1525 {
1526
1527   gStyle->SetOptTitle(0);
1528   gStyle->SetOptStat(0);
1529   gStyle->SetOptFit(0);
1530   CutSet = cutset;
1531   TString filename, simfilename;
1532   if(cutset==0){
1533     if(isPhos){
1534       if(isMC){
1535         filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
1536       }
1537       else{
1538         filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.PHOS.LHC10hPass2.Run139465.root";
1539       }
1540       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
1541     }
1542     else{
1543       //filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.LooseTrackMatchCuts.root";
1544       //simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.LooseTrackMatchCuts.root";
1545       if(isMC){
1546         filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
1547       }
1548       else{
1549         filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.root";
1550       }
1551       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
1552     }
1553   }
1554   if(cutset==1){
1555     if(isPhos){
1556       filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.PHOS350MeVCut.LHC10hPass2.Run139465.root";
1557       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS350MeVCut.LHC11a10a_bis.Run139465.root";
1558     }
1559     else{
1560       //filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.LooseTrackMatchCuts.root";
1561       //simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.LooseTrackMatchCuts.root";
1562       filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal350MeVCut.LHC10hPass2.Run139465.root";
1563       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal350MeVCut.LHC11a10a_bis.Run139465.root";
1564     }
1565   }
1566   if(cutset==2){
1567     if(isPhos){
1568       filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.PHOS400MeVCut.LHC10hPass2.Run139465.root";
1569       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS400MeVCut.LHC11a10a_bis.Run139465.root";
1570     }
1571     else{
1572       //filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal.LHC10hPass2.Run139465.LooseTrackMatchCuts.root";
1573       //simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.LooseTrackMatchCuts.root";
1574       filename = "rootFiles/LHC10hPass2/Et.ESD.realPbPb.EMCal400MeVCut.LHC10hPass2.Run139465.root";
1575       simfilename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal400MeVCut.LHC11a10a_bis.Run139465.root";
1576     }
1577   }
1578   cout<<"data file name = "<<filename<<endl;
1579   cout<<" sim file name = "<<simfilename<<endl;
1580   TString detector = "Emcal";
1581   Float_t scale = 360.0/40.0/1.2;//Azimuthal acceptance over eta range
1582   if(filename.Contains("PHOS")){
1583     detector = "Phos";
1584     //isPhos = kTRUE;
1585     scale = 360.0/60/0.24;//Azimuthal acceptance over eta range
1586   }
1587   TString det = detector;
1588   //gROOT->LoadMacro("macros/PlotSecondariesCorr.C");
1589   //PlotSecondariesCorr(simfilename,filename);
1590   ReadInData(filename,detector);
1591   ReadInNeutronCorrections();
1592   ReadInSecondaryCorrections();
1593   ReadInKaonCorrections();
1594   ReadMinEtCorrections();
1595   CalculateHadronCorrections(isPhos);
1596   ApplyCorrections(scale);
1597
1598   TCanvas *c1 = new TCanvas("c1","Corrections",500,400);
1599   c1->SetTopMargin(0.02);
1600   c1->SetRightMargin(0.02);
1601   c1->SetBorderSize(0);
1602   c1->SetFillColor(0);
1603   c1->SetFillColor(0);
1604   c1->SetBorderMode(0);
1605   c1->SetFrameFillColor(0);
1606   c1->SetFrameBorderMode(0);
1607   c1->SetLeftMargin(0.120968);
1608   //c1->SetRightMargin();
1609   //c1->SetTopMargin();
1610   c1->SetBottomMargin(0.15);
1611   c1->SetLogx();
1612   TGraphErrors *graphKaonCorrectionShort = GetKaonCorrectionGraphShort();
1613   graphKaonCorrectionShort->Draw("AP");
1614   //graphKaonCorrectionShort->SetMaximim(0.02);
1615   //graphKaonCorrectionShort->GetYaxis()->SetRange(1,graphKaonCorrectionShort->GetYaxis()->FindBin(0.02));
1616   float emcalmax = 0.006;
1617   graphKaonCorrectionShort->GetHistogram()->SetMaximum(emcalmax);
1618   //set scales the same within naive geometric scaling
1619   if(isPhos) graphKaonCorrectionShort->GetHistogram()->SetMaximum(emcalmax*0.24/1.2*60.0/40.0);
1620   graphKaonCorrectionShort->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
1621   graphKaonCorrectionShort->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
1622   graphKaonCorrectionShort->GetHistogram()->GetXaxis()->SetTitleSize(0.06);
1623   graphKaonCorrectionShort->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
1624   graphKaonCorrectionShort->GetXaxis()->SetTitle("N_{Ch}");
1625   graphKaonCorrectionShort->GetYaxis()->SetTitle("Correction/N_{Ch}");
1626   TGraphErrors *graphKaonCorrection = GetKaonCorrectionGraph();
1627   graphKaonCorrection->Draw("P same");
1628   TGraphErrors *graphSecondaryCorrectionShort = GetSecondaryCorrectionGraphShort();
1629   graphSecondaryCorrectionShort->Draw("P same");
1630   TGraphErrors *graphNeutronCorrectionShort = GetNeutronCorrectionGraphShort();
1631   graphNeutronCorrectionShort->Draw("P same");
1632   TGraphErrors *graphSecondaryCorrection = GetSecondaryCorrectionGraph();
1633   graphSecondaryCorrection->Draw("P same");
1634   TGraphErrors *graphNeutronCorrection = GetNeutronCorrectionGraph();
1635   graphNeutronCorrection->Draw("P same");
1636   TGraphErrors *graphHadronCorrection = GetHadronCorrectionGraph();
1637   graphHadronCorrection->Draw("P same");
1638   TGraphErrors *graphHadronCorrectionShort = GetHadronCorrectionGraphShort();
1639   graphHadronCorrectionShort->Draw("P same");
1640   TGraphErrors *graphTotalCorrection = GetTotalCorrectionPerNChGraph();
1641   graphTotalCorrection->Draw("P same");
1642   TLegend *leg = new TLegend(0.538306,0.756684,0.739919,0.97861);
1643   leg->SetFillStyle(0);
1644   leg->SetFillColor(0);
1645   leg->SetBorderSize(0);
1646   leg->SetTextSize(0.03);
1647   leg->SetTextSize(0.038682);
1648   leg->AddEntry(graphKaonCorrectionShort,"Kaon correction/N_{ch}","P");
1649   leg->AddEntry(graphNeutronCorrectionShort,"Neutron correction/N_{ch}","P");
1650   leg->AddEntry(graphSecondaryCorrectionShort,"Secondary correction/5.0 N_{ch}","P");
1651   leg->AddEntry(graphHadronCorrectionShort,"Hadron correction/5.0 N_{ch}","P");
1652   leg->AddEntry(graphTotalCorrection,"(E_{T}^{kaon}+E_{T}^{n}+E_{T}^{secondary}+E_{T}^{h})/20.0 N_{ch}","P");
1653   leg->Draw();
1654   TString corrplotname1 = "/tmp/CorrectionsVsNch"+detector+".eps";
1655   c1->SaveAs(corrplotname1.Data());
1656 //   TH1F *frame = new TH1F("frame","frame",1,0,2);
1657 //   frame->GetYaxis()->SetTitle("dE_{T}/d#eta");
1658 //   frame->GetXaxis()->SetTitle("N_{part}");
1659 //   //fPion->SetRange(0,2);
1660 //   frame->SetMinimum(0);
1661 //   frame->SetMaximum(10);
1662 //   TCanvas *c1a = new TCanvas("c1a","CorrectionsNoEffCorr",500,400);
1663 //   c1a->SetTopMargin(0.02);
1664 //   c1a->SetRightMargin(0.02);
1665 //   c1a->SetBorderSize(0);
1666 //   c1a->SetFillColor(0);
1667 //   c1a->SetFillColor(0);
1668 //   c1a->SetBorderMode(0);
1669 //   c1a->SetFrameFillColor(0);
1670 //   c1a->SetFrameBorderMode(0);
1671 //   c1a->SetLeftMargin(0.120968);
1672 //   //c1a->SetRightMargin();
1673 //   //c1a->SetTopMargin();
1674 //   c1a->SetBottomMargin(0.15);
1675 //   c1a->SetLogx();
1676 //   TGraphErrors *graphKaonCorrectionShortNoEffCorr = GetKaonCorrectionGraphShortNoEffCorr();
1677 //   graphKaonCorrectionShortNoEffCorr->Draw("AP");
1678 //   //graphKaonCorrectionShort->SetMaximim(0.02);
1679 //   //graphKaonCorrectionShort->GetYaxis()->SetRange(1,graphKaonCorrectionShort->GetYaxis()->FindBin(0.02));
1680 //   float emcalmax = 0.006;
1681 //   graphKaonCorrectionShortNoEffCorr->GetHistogram()->SetMaximum(emcalmax);
1682 //   //set scales the same within naive geometric scaling
1683 //   if(isPhos) graphKaonCorrectionShortNoEffCorr->GetHistogram()->SetMaximum(emcalmax*0.24/1.2*60.0/40.0);
1684 //   graphKaonCorrectionShortNoEffCorr->GetHistogram()->GetXaxis()->SetLabelSize(0.04);
1685 //   graphKaonCorrectionShortNoEffCorr->GetHistogram()->GetYaxis()->SetLabelSize(0.04);
1686 //   graphKaonCorrectionShortNoEffCorr->GetHistogram()->GetXaxis()->SetTitleSize(0.06);
1687 //   graphKaonCorrectionShortNoEffCorr->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
1688 //   graphKaonCorrectionShortNoEffCorr->GetXaxis()->SetTitle("N_{Ch}");
1689 //   graphKaonCorrectionShortNoEffCorr->GetYaxis()->SetTitle("Correction/N_{Ch}");
1690 //   TGraphErrors *graphKaonCorrectionNoEffCorr = GetKaonCorrectionGraphNoEffCorr();
1691 //   graphKaonCorrectionNoEffCorr->Draw("P same");
1692 //   TGraphErrors *graphSecondaryCorrectionShortNoEffCorr = GetSecondaryCorrectionGraphShortNoEffCorr();
1693 //   graphSecondaryCorrectionShortNoEffCorr->Draw("P same");
1694 //   TGraphErrors *graphNeutronCorrectionShortNoEffCorr = GetNeutronCorrectionGraphShortNoEffCorr();
1695 //   graphNeutronCorrectionShortNoEffCorr->Draw("P same");
1696 //   TGraphErrors *graphSecondaryCorrectionNoEffCorr = GetSecondaryCorrectionGraphNoEffCorr();
1697 //   graphSecondaryCorrectionNoEffCorr->Draw("P same");
1698 //   TGraphErrors *graphNeutronCorrectionNoEffCorr = GetNeutronCorrectionGraphNoEffCorr();
1699 //   graphNeutronCorrectionNoEffCorr->Draw("P same");
1700 //   TGraphErrors *graphHadronCorrectionNoEffCorr = GetHadronCorrectionGraphNoEffCorr();
1701 //   graphHadronCorrectionNoEffCorr->Draw("P same");
1702 //   TGraphErrors *graphHadronCorrectionShortNoEffCorr = GetHadronCorrectionGraphShortNoEffCorr();
1703 //   graphHadronCorrectionShortNoEffCorr->Draw("P same");
1704 //   TGraphErrors *graphTotalCorrectionNoEffCorr = GetTotalCorrectionPerNChGraphNoEffCorr();
1705 //   graphTotalCorrectionNoEffCorr->Draw("P same");
1706 //   TLegend *leg2 = new TLegend(0.538306,0.756684,0.739919,0.97861);
1707 //   leg2->SetFillStyle(0);
1708 //   leg2->SetFillColor(0);
1709 //   leg2->SetBorderSize(0);
1710 //   leg2->SetTextSize(0.03);
1711 //   leg2->SetTextSize(0.038682);
1712 //   leg2->AddEntry(graphKaonCorrectionShortNoEffCorr,"Kaon correction/N_{ch}","P");
1713 //   leg2->AddEntry(graphNeutronCorrectionShortNoEffCorr,"Neutron correction/N_{ch}","P");
1714 //   leg2->AddEntry(graphSecondaryCorrectionShortNoEffCorr,"Secondary correction/5.0 N_{ch}","P");
1715 //   leg2->AddEntry(graphHadronCorrectionShortNoEffCorr,"Hadron correction/5.0 N_{ch}","P");
1716 //   leg2->AddEntry(graphTotalCorrectionNoEffCorr,"(E_{T}^{kaon}+E_{T}^{n}+E_{T}^{secondary}+E_{T}^{h})/20.0 N_{ch}","P");
1717 //   leg2->Draw();
1718 //   TString corrplotname1 = "/tmp/CorrectionsVsNch"+detector+"NoEffCorr.eps";
1719 //   c1a->SaveAs(corrplotname1.Data());
1720   TCanvas *c2 = new TCanvas("c2","Min Et Correction",600,400);
1721   c2->SetTopMargin(0.02);
1722   c2->SetRightMargin(0.02);
1723   c2->SetBorderSize(0);
1724   c2->SetFillColor(0);
1725   c2->SetFillColor(0);
1726   c2->SetBorderMode(0);
1727   c2->SetFrameFillColor(0);
1728   c2->SetFrameBorderMode(0);
1729   c2->SetLogx();
1730   TGraphErrors *graphMinEtCorrectionShort = GetMinEtCorrectionGraphShort();
1731   graphMinEtCorrectionShort->GetXaxis()->SetTitle("f_{ETmin}");
1732   graphMinEtCorrectionShort->GetYaxis()->SetTitle("Correction/N_{Ch}");
1733 //   graphMinEtCorrectionShort->GetHistogram()->SetMaximum(1.0);
1734 //   graphMinEtCorrectionShort->GetHistogram()->SetMinimum(0.0);
1735   graphMinEtCorrectionShort->Draw("AP");
1736   graphKaonCorrectionShort->GetXaxis()->SetTitle("N_{Ch}");
1737   graphKaonCorrectionShort->GetYaxis()->SetTitle("f_{EtMin}");
1738   TGraphErrors *graphMinEtCorrection = GetMinEtCorrectionGraph();
1739   graphMinEtCorrection->Draw("P same");
1740   TString corrplotname2 = "/tmp/EtMinVsNch"+detector+".eps";
1741   c2->SaveAs(corrplotname2.Data());
1742
1743 //   TCanvas *c3 = new TCanvas("c3","Kaon Correction",600,400);
1744 //   c3->SetTopMargin(0.02);
1745 //   c3->SetRightMargin(0.02);
1746 //   c3->SetBorderSize(0);
1747 //   c3->SetFillColor(0);
1748 //   c3->SetFillColor(0);
1749 //   c3->SetBorderMode(0);
1750 //   c3->SetFrameFillColor(0);
1751 //   c3->SetFrameBorderMode(0);
1752 //   c3->SetLogx();
1753
1754 //   TGraphErrors *graphKaonGraph = GetKaonGraph();
1755 //   graphKaonGraph->Draw("AP");
1756 //   //graphKaonGraph->GetHistogram()->SetMaximum();
1757 //   graphKaonGraph->GetHistogram()->SetMinimum(0.0);
1758 //   graphKaonGraph->GetXaxis()->SetTitle("N_{Ch}");
1759 //   graphKaonGraph->GetYaxis()->SetTitle("val/N_{Ch}");
1760 //   graphKaonCorrectionShort->Draw("P same");
1761 //   TGraphErrors *graphKaonEtGraph = GetKaonEtGraph();
1762 //   graphKaonEtGraph->Draw("P same");
1763 //   TLegend *leg3 = new TLegend(0.602349,0.362903,0.805369,0.586022);
1764 //   leg3->SetFillStyle(0);
1765 //   leg3->SetFillColor(0);
1766 //   leg3->SetBorderSize(0);
1767 //   leg3->SetTextSize(0.03);
1768 //   leg3->SetTextSize(0.038682);
1769 //   leg3->AddEntry(graphKaonGraph,"N_{K^{#pm}}/N_{Ch}","P");
1770 //   leg3->AddEntry(graphKaonCorrectionShort,"E_{T}^{kaon}","P");
1771 //   leg3->AddEntry(graphKaonEtGraph,"E_{T}^{K^{#pm}}/N_{Ch}","P");
1772 //   leg3->Draw();
1773 //   TString corrplotname3 = "/tmp/KaonVsNch"+detector+".eps";
1774 //   c3->SaveAs(corrplotname3.Data());
1775
1776
1777   TCanvas *c4 = new TCanvas("c4", "dE_{T}/d#eta#frac{1}{0.5*N_{part}} [GeV]",700, 600);
1778   c4->SetTopMargin(0.02);
1779   c4->SetRightMargin(0.02);
1780   c4->SetBorderSize(0);
1781   c4->SetFillColor(0);
1782   c4->SetFillColor(0);
1783   c4->SetBorderMode(0);
1784   c4->SetFrameFillColor(0);
1785   c4->SetFrameBorderMode(0);
1786   TGraphErrors *graphEt = GetEtGraph();
1787   graphEt->GetHistogram()->GetXaxis()->SetTitle("N_{part}");
1788   graphEt->GetHistogram()->GetYaxis()->SetTitle("(E_{T}/(N_{part}/2)");
1789   graphEt->GetHistogram()->SetMinimum(0.0);
1790   graphEt->GetHistogram()->SetMaximum(3.0);
1791   graphEt->Draw("AP");
1792   TGraphErrors *graphPartialCorr = GetPartialCorrEtPerNPartPairGraph();
1793   //graphPartialCorr->Draw("P same");
1794   TGraphErrors *graphPionEt = GetPionEtGraph();
1795   graphPionEt->Draw("P same");
1796   TGraphErrors *graphPionEmEt = GetPionEmEtGraph();
1797   graphPionEmEt->Draw("P same");
1798   TGraphErrors *graphEtFormulaC = GetEtGraphFormulaC();
1799   //graphEtFormulaC->Draw("P same");
1800   TGraphErrors *graphEtFormulaB = GetEtGraphFormulaB();
1801   //graphEtFormulaB->Draw("P same");
1802   TLegend *leg2 = new TLegend(0.607383,0.704301,0.810403,0.927419);
1803   leg2->SetFillStyle(0);
1804   leg2->SetFillColor(0);
1805   leg2->SetBorderSize(0);
1806   leg2->SetTextSize(0.03);
1807   leg2->SetTextSize(0.038682);
1808   leg2->AddEntry(graphEt,"E_{T}^{EM}","P");
1809   leg2->AddEntry(graphPionEt,"E_{T}^{#pi^{#pm}}","P");
1810   leg2->AddEntry(graphPionEmEt,"E_{T}^{EM calc. from #pi}","P");
1811   //leg2->AddEntry(graphEtFormulaC,"E_{T}^{EM} Formula C","P");
1812   //leg2->AddEntry(graphEtFormulaB,"E_{T}^{EM} Formula B","P");
1813   leg2->Draw();
1814   TString corrplotname4 = "/tmp/EtVsNpart"+detector+".eps";
1815   c4->SaveAs(corrplotname4.Data());
1816
1817 //   TCanvas *c5 = new TCanvas("c5", "Raw Et Vs NCh",1200, 400);
1818 //   c5->SetTopMargin(0.02);
1819 //   c5->SetRightMargin(0.02);
1820 //   c5->SetBorderSize(0);
1821 //   c5->SetFillColor(0);
1822 //   c5->SetFillColor(0);
1823 //   c5->SetBorderMode(0);
1824 //   c5->SetFrameFillColor(0);
1825 //   c5->SetFrameBorderMode(0);
1826 //   c5->Divide(2);
1827 //   c5->cd(1);
1828 //   TGraphErrors *partialCorrEtPerNCh = GetPartialCorrEtPerNChGraph();
1829 //   partialCorrEtPerNCh->GetHistogram()->SetMaximum(0.1*7.5/scale);//so PHOS and EMCal scales are similar...  partialCorrEtPerNCh->GetHistogram()->SetMinimum(0.0);
1830 //   partialCorrEtPerNCh->GetHistogram()->GetXaxis()->SetTitle("N_{Ch}");
1831 //   partialCorrEtPerNCh->GetHistogram()->GetYaxis()->SetTitle("(E_{T}/#epsilon f_{nonlin})/N_{Ch}");
1832 //   partialCorrEtPerNCh->Draw("AP");
1833 //   TGraphErrors *graphRawNoEffCorrCh = GetRawEtNoEffCorrPerNChGraph();
1834 //   graphRawNoEffCorrCh->Draw("P same");
1835 //   TGraphErrors *graphRawAllNoEffCorrCh = GetRawEtAllNoEffCorrPerNChGraph();
1836 //   graphRawAllNoEffCorrCh->Draw("P same");
1837 //   TGraphErrors *graphRawAllCh = GetRawEtAllPerNChGraph();
1838 //   graphRawAllCh->Draw("P same");
1839 //   TLegend *leg5 = new TLegend(0.416667,0.173077,0.62069,0.396853);
1840 //   leg5->SetFillStyle(0);
1841 //   leg5->SetFillColor(0);
1842 //   leg5->SetBorderSize(0);
1843 //   leg5->SetTextSize(0.03);
1844 //   leg5->SetTextSize(0.038682);
1845 //   leg5->AddEntry(partialCorrEtPerNCh,"E_{T}#delta_{TM}/#epsilon f_{nonlin}","P");
1846 //   leg5->AddEntry(graphRawNoEffCorrCh,"E_{T}#delta_{TM}/ f_{nonlin}","P");
1847 //   leg5->AddEntry(graphRawAllCh,"E_{T}/#epsilon f_{nonlin}","P");
1848 //   leg5->AddEntry(graphRawAllNoEffCorrCh,"E_{T}/f_{nonlin}","P");
1849 //   leg5->Draw();
1850 //   TString corrplotname5 = "/tmp/RawEtVsNch"+detector+".eps";
1851 //   c5->SaveAs(corrplotname5.Data());
1852
1853 // //   TCanvas *c6 = new TCanvas("c6", "Raw Et Vs NCl",700, 600);
1854 // //   c6->SetTopMargin(0.02);
1855 // //   c6->SetRightMargin(0.02);
1856 // //   c6->SetBorderSize(0);
1857 // //   c6->SetFillColor(0);
1858 // //   c6->SetFillColor(0);
1859 // //   c6->SetBorderMode(0);
1860 // //   c6->SetFrameFillColor(0);
1861 // //   c6->SetFrameBorderMode(0);
1862
1863 //   c5->cd(2);
1864 //   TGraphErrors *graphRawAllCl = GetRawEtAllPerNClGraph();
1865 //   graphRawAllCl->GetHistogram()->GetXaxis()->SetTitle("N_{Ch}");
1866 //   graphRawAllCl->GetHistogram()->GetYaxis()->SetTitle("(E_{T}/#epsilon f_{nonlin})/N_{Ch}");
1867 //   graphRawAllCl->Draw("AP");
1868 //   TGraphErrors *partialCorrEtPerNCl = GetPartialCorrEtPerNClGraph();
1869 //   partialCorrEtPerNCl->GetHistogram()->GetXaxis()->SetTitle("N_{Ch}");
1870 //   partialCorrEtPerNCl->GetHistogram()->GetYaxis()->SetTitle("(E_{T}/#epsilon f_{nonlin})/N_{Ch}");
1871 //   partialCorrEtPerNCl->Draw("P same");
1872 //   TGraphErrors *graphRawNoEffCorrCl = GetRawEtNoEffCorrPerNClGraph();
1873 //   graphRawNoEffCorrCl->Draw("P same");
1874 //   TGraphErrors *graphRawAllNoEffCorrCl = GetRawEtAllNoEffCorrPerNClGraph();
1875 //   graphRawAllNoEffCorrCl->Draw("P same");
1876 //   TLegend *leg6 = new TLegend(0.416667,0.173077,0.62069,0.396853);
1877 //   leg6->SetFillStyle(0);
1878 //   leg6->SetFillColor(0);
1879 //   leg6->SetBorderSize(0);
1880 //   leg6->SetTextSize(0.03);
1881 //   leg6->SetTextSize(0.038682);
1882 //   leg6->AddEntry(partialCorrEtPerNCl,"E_{T}#delta_{TM}/#epsilon f_{nonlin}","P");
1883 //   leg6->AddEntry(graphRawNoEffCorrCl,"E_{T}#delta_{TM}/ f_{nonlin}","P");
1884 //   leg6->AddEntry(graphRawAllCl,"E_{T}/#epsilon f_{nonlin}","P");
1885 //   leg6->AddEntry(graphRawAllNoEffCorrCl,"E_{T}/f_{nonlin}","P");
1886 //   leg6->Draw();
1887 //   TString corrplotname5 = "/tmp/RawEtVsNcl"+detector+".eps";
1888 //   c5->SaveAs(corrplotname5.Data());
1889
1890 //   TCanvas *c7 = new TCanvas("c7", "Matched tracks",1200, 400);
1891 //   c7->SetTopMargin(0.02);
1892 //   c7->SetRightMargin(0.02);
1893 //   c7->SetBorderSize(0);
1894 //   c7->SetFillColor(0);
1895 //   c7->SetFillColor(0);
1896 //   c7->SetBorderMode(0);
1897 //   c7->SetFrameFillColor(0);
1898 //   c7->SetFrameBorderMode(0);
1899 //   c7->Divide(2);
1900 //   c7->cd(1);
1901 //   TGraphErrors *graphTotalPerNCh = GetTotalTracksPerNChGraph();
1902 //   graphTotalPerNCh->GetHistogram()->GetXaxis()->SetTitle("N_{Ch}");
1903 //   graphTotalPerNCh->GetHistogram()->GetYaxis()->SetTitle("Number of tracks/N_{Ch}");
1904 //   graphTotalPerNCh->GetHistogram()->SetMaximum(0.06*7.5/scale);//so PHOS and EMCal scales are similar...
1905 //   graphTotalPerNCh->Draw("AP");
1906 //   TGraphErrors *graphMatchedPerNCh = GetMatchedTracksPerNChGraph();
1907 //   graphMatchedPerNCh->Draw("P same");
1908 //   TGraphErrors *graphNotMatchedPerNCh = GetNotMatchedTracksPerNChGraph();
1909 //   graphNotMatchedPerNCh->Draw("P same");
1910 //   TLegend *leg7 = new TLegend(0.176003,0.637152,0.379808,0.86208);
1911 //   leg7->SetFillStyle(0);
1912 //   leg7->SetFillColor(0);
1913 //   leg7->SetBorderSize(0);
1914 //   leg7->SetTextSize(0.03);
1915 //   leg7->SetTextSize(0.038682);
1916 //   leg7->AddEntry(graphMatchedPerNCh,"Number of tracks matched to clusters","P");
1917 //   leg7->AddEntry(graphNotMatchedPerNCh,"Number of tracks not matched to clusters","P");
1918 //   leg7->AddEntry(graphTotalPerNCh,"Total number of tracks","P");
1919 //   leg7->Draw();
1920 //   c7->cd(2);
1921 //   TGraphErrors *graphTotalPerNCl = GetTotalTracksPerNClGraph();
1922 //   graphTotalPerNCl->GetHistogram()->GetXaxis()->SetTitle("N_{Cl}");
1923 //   graphTotalPerNCl->GetHistogram()->GetYaxis()->SetTitle("Number of tracks/N_{Cl}");
1924 //   graphTotalPerNCl->GetHistogram()->SetMaximum(0.65);
1925 //   graphTotalPerNCl->Draw("AP");
1926 //   TGraphErrors *graphMatchedPerNCl = GetMatchedTracksPerNClGraph();
1927 //   graphMatchedPerNCl->Draw("P same");
1928 //   TGraphErrors *graphNotMatchedPerNCl = GetNotMatchedTracksPerNClGraph();
1929 //   graphNotMatchedPerNCl->Draw("P same");
1930 //   TLine *line = new TLine(graphTotalPerNCl->GetHistogram()->GetXaxis()->GetBinLowEdge(1),0.5,graphTotalPerNCl->GetHistogram()->GetXaxis()->GetBinLowEdge(graphTotalPerNCl->GetHistogram()->GetXaxis()->GetNbins()),0.5);
1931 //   line->Draw();
1932 //   TString corrplotname7 = "/tmp/TrackMatchCrossCheck"+detector+".eps";
1933 //   c7->SaveAs(corrplotname7.Data());
1934
1935   TCanvas *c8 = new TCanvas("c8","Signal Fraction",600,400);
1936   c8->SetTopMargin(0.02);
1937   c8->SetRightMargin(0.02);
1938   c8->SetBorderSize(0);
1939   c8->SetFillColor(0);
1940   c8->SetFillColor(0);
1941   c8->SetBorderMode(0);
1942   c8->SetFrameFillColor(0);
1943   c8->SetFrameBorderMode(0);
1944   //c8->SetLogx();
1945   TGraphErrors *graphSignalFraction = GetSignalFractionGraph();
1946   graphSignalFraction->GetXaxis()->SetTitle("f_{ETmin}");
1947   graphSignalFraction->GetYaxis()->SetTitle("Correction/N_{Ch}");
1948   graphSignalFraction->GetHistogram()->SetMaximum(0.4);//so PHOS and EMCal scales are similar...
1949   graphSignalFraction->GetHistogram()->SetMinimum(0.0);//so PHOS and EMCal scales are similar...
1950 //   graphSignalFraction->GetHistogram()->SetMaximum(1.0);
1951 //   graphSignalFraction->GetHistogram()->SetMinimum(0.0);
1952   graphSignalFraction->Draw("AP");
1953   TString corrplotname8 = "/tmp/SignalFraction"+detector+".eps";
1954   c8->SaveAs(corrplotname8.Data());
1955
1956   TCanvas *c9 = new TCanvas("c9","Signal Fraction",600,400);
1957   c9->SetTopMargin(0.02);
1958   c9->SetRightMargin(0.02);
1959   c9->SetBorderSize(0);
1960   c9->SetFillColor(0);
1961   c9->SetFillColor(0);
1962   c9->SetBorderMode(0);
1963   c9->SetFrameFillColor(0);
1964   c9->SetFrameBorderMode(0);
1965   c9->SetLogy();
1966   c9->SetLogx();
1967   int nbins = 15;
1968   if(isPhos) nbins = 8;
1969   //nbins = 20;
1970 //   for(bin = 1; bin<=nbins;bin++){
1971 //     dataFits[bin] = new TF1(Form("Fit%i",bin),"[0]*exp(-(x-[1])*(x-[1])/[2])",0,250);
1972 //     ((TF1*)dataFits[bin])->SetParameter(0,1e3);
1973 //     ((TF1*)dataFits[bin])->SetParameter(1,partialCorrEtValues[bin-1]);
1974 //     ((TF1*)dataFits[bin])->SetParameter(2,20);
1975 //     ((TH1D*)rawEt[bin])->Fit(((TF1*)dataFits[bin]));
1976 //   }
1977   int bin = 1;
1978   ((TH1D*)rawEt[bin])->GetXaxis()->SetTitle("E_{T}^{raw}");
1979   if(isPhos){((TH1D*)rawEt[bin])->GetXaxis()->SetRange(1,((TH1D*)rawEt[bin])->GetXaxis()->FindBin(100));}
1980   else{((TH1D*)rawEt[bin])->GetXaxis()->SetRange(1,((TH1D*)rawEt[bin])->GetXaxis()->FindBin(200));}
1981   ((TH1D*)rawEt[bin])->Draw();
1982   ((TH1D*)rawEt[bin])->SetMaximum(5e4);
1983   for(bin = 1; bin<=nbins;bin+=1){
1984     ((TH1D*)rawEt[bin])->Draw("same");
1985     ((TH1D*)rawEt[bin])->SetLineColor(colors[bin-1]);
1986     ((TH1D*)rawEt[bin])->SetMarkerColor(colors[bin-1]);
1987     //((TF1*)dataFits[bin])->SetLineColor(colors[bin-1]);
1988     //Float_t mostlikely = GetMostLikelyValue((TH1D*)rawEt[bin]);
1989     //cout<<"cb "<<bin-1<<" average: "<<partialCorrEtValues[bin-1]<<" most likely "<<mostlikely<<" fit "<<((TF1*)dataFits[bin])->GetParameter(1)<<endl;
1990   }
1991   TString name9 = "/tmp/ETDistribution"+detector+".png";
1992   c9->SaveAs(name9.Data());
1993   //Create correction file with fractional corrections so that I can use them for cross checks
1994   ofstream myfile;
1995   TString texfilename = "SignalFrac"+detector+".dat";
1996   myfile.open (texfilename.Data());
1997   //neutron, hadron, kaon, secondaries
1998   for(int i=0;i<20;i++){
1999     myfile <<signalFraction[i]<<"\t";
2000     myfile <<signalFractionError[i]<<"\t";
2001     myfile <<endl;
2002   }
2003   myfile.close();
2004
2005   //Create correction file with all corrections so that I can use them for cross checks
2006   ofstream myfile;
2007   TString texfilename = "allcorrections"+detector+".dat";
2008   myfile.open (texfilename.Data());
2009   //neutron, hadron, kaon, secondaries
2010   for(int i=0;i<20;i++){
2011     myfile <<neutronCorr[i]<<"\t";
2012     myfile <<hadCorr[i]<<"\t";
2013     myfile <<kaonCorr[i]<<"\t";
2014     myfile <<secondaryCorr[i]<<endl;
2015   }
2016   myfile.close();
2017   //from spectra calc
2018   //0.232 $\pm$ 0.029
2019   Float_t fem = 0.232;
2020   Float_t femerr = 0.029;
2021   ofstream myfile2;
2022   ofstream myfile3;
2023   //EmEtFromSpectra.dat
2024   TString texfilename = "EmEtFrom"+detector+".dat";
2025   TString texfilename3 = "TotEtFrom"+detector+".dat";
2026   myfile2.open (texfilename.Data());
2027   myfile3.open (texfilename3.Data());
2028   Int_t maxcb = 10;
2029   if(isPhos) maxcb = 6;
2030   //neutron, hadron, kaon, secondaries
2031   for(int i=0;i<maxcb;i++){
2032     int cb1 = i*5;
2033     int cb2 = (i+1)*5;
2034     myfile2 <<cb1<<"\t"<<cb2<<"\t";
2035     myfile2 <<corrEtValues[i]<<"\t";
2036     myfile2 <<corrEtError[i]<<"\t";
2037     myfile2 <<endl;
2038
2039     Float_t value = corrEtValues[i]/fem;
2040     Float_t error = value * TMath::Sqrt(TMath::Power(corrEtError[i]/corrEtValues[i],2)+TMath::Power(femerr/fem,2));
2041     myfile3 <<cb1<<"\t"<<cb2<<"\t";
2042     myfile3 <<value<<"\t";
2043     myfile3 <<error<<"\t";
2044     myfile3 <<endl;
2045   }
2046   myfile2.close();
2047   myfile3.close();
2048
2049   WriteLatex();
2050 }
2051
2052
2053 void WriteLatex(){
2054   TString detector = "Emcal";
2055     string inline;
2056     //FinalemetsEmcal.dat
2057     TString finalemetInfileName = "EmEtFrom"+detector+".dat";
2058     ifstream myfinalemetfile3 (finalemetInfileName.Data());
2059     Float_t value = 0;
2060     Float_t error = 0;
2061     Float_t junk1 = 0;
2062     Float_t junk2 = 0;
2063     Int_t i=0;
2064     if (myfinalemetfile3.is_open()){
2065       while ( myfinalemetfile3.good() )
2066         {
2067           getline (myfinalemetfile3,inline);
2068           istringstream tmp(inline);
2069           tmp >> junk1;
2070           tmp >> junk2;
2071           tmp >> value;
2072           tmp >> error;
2073           if(i<20){
2074             //cout<<"value "<<value<<" +/- "<<error<<endl;
2075             finalemetCorrEmcal[i] = value;
2076             finalemetErrorEmcal[i] = error;
2077           }
2078           i++;
2079         }
2080         myfinalemetfile3.close();
2081     }
2082
2083     detector = "Phos";
2084     finalemetInfileName = "EmEtFrom"+detector+".dat";
2085     ifstream myfinalemetfile4 (finalemetInfileName.Data());
2086     Float_t value = 0;
2087     Float_t error = 0;
2088     Int_t i=0;
2089     if (myfinalemetfile4.is_open()){
2090       while ( myfinalemetfile4.good() )
2091         {
2092           getline (myfinalemetfile4,inline);
2093           istringstream tmp(inline);
2094           tmp >> junk1;
2095           tmp >> junk2;
2096           tmp >> value;
2097           tmp >> error;
2098           if(i<20){
2099             finalemetCorrPhos[i] = value;
2100             finalemetErrorPhos[i] = error;
2101           }
2102           i++;
2103         }
2104         myfinalemetfile4.close();
2105     }
2106     for(int i=0;i<20;i++){
2107       TString line = Form("%i-%i & %2.3f $\\pm$ %2.3f & %2.3f $\\pm$ %2.3f \\\\",i*5,(i+1)*5,finalemetCorrPhos[i],finalemetErrorPhos[i],finalemetCorrEmcal[i],finalemetErrorEmcal[i]);
2108       cout<<line.Data()<<endl;
2109     }
2110
2111
2112
2113
2114
2115
2116
2117   detector = "Emcal";
2118     string inline;
2119     //FinaltotaletsEmcal.dat
2120     TString finaltotaletInfileName = "TotEtFrom"+detector+".dat";
2121     ifstream myfinaltotaletfile5 (finaltotaletInfileName.Data());
2122      value = 0;
2123      error = 0;
2124      junk1 = 0;
2125      junk2 = 0;
2126     Int_t i=0;
2127     if (myfinaltotaletfile5.is_open()){
2128       while ( myfinaltotaletfile5.good() )
2129         {
2130           getline (myfinaltotaletfile5,inline);
2131           istringstream tmp(inline);
2132           tmp >> junk1;
2133           tmp >> junk2;
2134           tmp >> value;
2135           tmp >> error;
2136           if(i<20){
2137             //cout<<"value "<<value<<" +/- "<<error<<endl;
2138             finaltotaletCorrEmcal[i] = value;
2139             finaltotaletErrorEmcal[i] = error;
2140           }
2141           i++;
2142         }
2143         myfinaltotaletfile5.close();
2144     }
2145
2146     detector = "Phos";
2147     finaltotaletInfileName = "TotEtFrom"+detector+".dat";
2148     ifstream myfinaltotaletfile6 (finaltotaletInfileName.Data());
2149      value = 0;
2150      error = 0;
2151     Int_t i=0;
2152     if (myfinaltotaletfile6.is_open()){
2153       while ( myfinaltotaletfile6.good() )
2154         {
2155           getline (myfinaltotaletfile6,inline);
2156           istringstream tmp(inline);
2157           tmp >> junk1;
2158           tmp >> junk2;
2159           tmp >> value;
2160           tmp >> error;
2161           if(i<20){
2162             finaltotaletCorrPhos[i] = value;
2163             finaltotaletErrorPhos[i] = error;
2164           }
2165           i++;
2166         }
2167         myfinaltotaletfile6.close();
2168     }
2169     for(int i=0;i<20;i++){
2170       TString line = Form("%i-%i & %2.3f $\\pm$ %2.3f & %2.3f $\\pm$ %2.3f \\\\",i*5,(i+1)*5,finaltotaletCorrPhos[i],finaltotaletErrorPhos[i],finaltotaletCorrEmcal[i],finaltotaletErrorEmcal[i]);
2171       cout<<line.Data()<<endl;
2172     }
2173
2174
2175
2176
2177
2178
2179 }