Peter Parfenov лет назад: 4
Родитель
Сommit
e75af1c70a
1 измененных файлов с 58 добавлено и 23 удалено
  1. 58 23
      FitCent.cpp

+ 58 - 23
FitCent.cpp

@@ -55,8 +55,8 @@ Double_t ftEtaN(Double_t *x, Double_t *par)
     Double_t xx = x[0];
     // Function
     Double_t mean_n = ftMeanN(&x[0], &par[1]);
-    Double_t sigma  = ftSigmN(&x[0], &par[0]);
-    Double_t func = 2 * TMath::Power(1+TMath::Erf(mean_n/(TMath::Sqrt(2)*sigma)), -1);
+    Double_t sigma = ftSigmN(&x[0], &par[0]);
+    Double_t func = 2 * TMath::Power(1 + TMath::Erf(mean_n / (TMath::Sqrt(2) * sigma)), -1);
 
     return func;
 }
@@ -95,7 +95,7 @@ Double_t ftProbIntegrandEta(Double_t *x, Double_t *par)
     // Function
     Double_t mean_n = ftMeanN(&x[0], &par[2]);
     Double_t sigma = ftSigmN(&x[0], &par[1]);
-    Double_t eta   = ftEtaN(&x[0], &par[1]);
+    Double_t eta = ftEtaN(&x[0], &par[1]);
     Double_t func = eta / (sigma * TMath::Sqrt(2 * TMath::Pi())) *
                     TMath::Exp(-TMath::Power((n - mean_n), 2) / (2 * sigma * sigma));
 
@@ -189,9 +189,9 @@ Double_t ftProbCb(Double_t *x, Double_t *par)
     // Function
     TF1 *f = new TF1("f", ftProbIntegrand, 0, 1, 6);
     f->SetParameters(n, sig0, n0, a1, a2, a3);
-    Double_t Denom = f->Integral(0,1);
-    Double_t Nom   = f->Eval(xx);
-    Double_t func  = Nom/Denom;
+    Double_t Denom = f->Integral(0, 1);
+    Double_t Nom = f->Eval(xx);
+    Double_t func = Nom / Denom;
 
     return func;
 }
@@ -207,12 +207,12 @@ Double_t ftProbB(Double_t *x, Double_t *par)
     Double_t a2 = par[5];
     Double_t a3 = par[6];
     // Variables
-    Double_t b  = x[0];
-    Double_t cb = TMath::Pi()*TMath::Power(b,2)/sigInel;
+    Double_t b = x[0];
+    Double_t cb = TMath::Pi() * TMath::Power(b, 2) / sigInel;
     // Function
     TF1 *f = new TF1("f", ftProbCb, 0, 1, 6);
     f->SetParameters(n, sig0, n0, a1, a2, a3);
-    Double_t func = 2*TMath::Pi()*b/sigInel*f->Eval(cb);
+    Double_t func = 2 * TMath::Pi() * b / sigInel * f->Eval(cb);
 
     return func;
 }
@@ -232,6 +232,13 @@ int main(int argc, char **argv)
     // Number of iterations during fitting procedure
     int NumOfIterations = 5;
 
+    // Reference statistics for hBimp
+    int RefStat = 1e6;
+
+    const TString centName[10] = {"0-10%", "10-20%", "20-30%", "30-40%",
+                                  "40-50%", "50-60%", "60-70%", "70-80%",
+                                  "80-90%", "90-100%"};
+
     // Parsing arguments
     TString inputfile, outputfile;
     Bool_t isRoot = false;
@@ -366,20 +373,43 @@ int main(int argc, char **argv)
     // Calculate probability functions for 10% step centrality
     TF1 *probFunc[10];
     TF1 *probFuncB[10];
-    Double_t n_test;
-    const Double_t sigma_Inelastic = 777.; // [fm^2] for Pb+Pb
-    for (int icent=0; icent<10; icent++)
+    TF1 *funcTmp;
+    TH1F *hBimp[10];
+    TH1F *hBimpTotal = new TH1F(Form("hBimpTotal"), Form("%s","0-100%"), 200, 0, 20);
+    Double_t n_test, n1, n2;
+    const Double_t sigma_Inelastic = 685.; // [fm^2] for Au+Au
+    // const Double_t sigma_Inelastic = 777.; // [fm^2] for Pb+Pb
+    for (int icent = 0; icent < 10; icent++)
     {
-        probFunc[icent] = new TF1(Form("probFunc_cent%i",icent), ftProbCb, 0, 1, 6);
-        probFuncB[icent] = new TF1(Form("probFuncB_cent%i",icent), ftProbB, 0, 20, 7);
-        n_test = centFunc->GetX(icent/10.+0.05);
-        probFunc[icent]->SetParameter(0,n_test);
-        probFuncB[icent]->SetParameter(0,sigma_Inelastic);
-        probFuncB[icent]->SetParameter(1,n_test);
-        for (int ipar=0; ipar<5; ipar++)
+        // std::cout << "CHECKIN' " << icent << std::endl;
+        probFunc[icent] = new TF1(Form("probFunc_cent%i", icent), ftProbCb, 0, 1, 6);
+        probFuncB[icent] = new TF1(Form("probFuncB_cent%i", icent), ftProbB, 0, 20, 7);
+        hBimp[icent] = new TH1F(Form("hBimp_cent%i", icent), Form("%s", centName[icent].Data()), 200, 0, 20);
+        for (int jcent = 0; jcent < 10; jcent++)
         {
-            probFunc[icent]->SetParameter(ipar+1,fitFunc->GetParameter(ipar));
-            probFuncB[icent]->SetParameter(ipar+2,fitFunc->GetParameter(ipar));
+            // std::cout << "\tCentrality = " << icent*10+jcent << std::endl;
+            if ((icent + jcent) == 0 || (icent * 10 + jcent) >= 98)
+                continue;
+            n_test = centFunc->GetX(icent / 10. + jcent / 100.);
+            funcTmp = new TF1(Form("funcTmp"), ftProbB, 0, 20, 7);
+            funcTmp->SetParameter(0, sigma_Inelastic);
+            funcTmp->SetParameter(1, n_test);
+            for (int ipar = 0; ipar < 5; ipar++)
+            {
+                funcTmp->SetParameter(ipar + 2, fitFunc->GetParameter(ipar));
+            }
+            hBimp[icent]->FillRandom("funcTmp", RefStat);
+            hBimpTotal->FillRandom("funcTmp", RefStat);
+            delete funcTmp;
+        }
+        n_test = centFunc->GetX(icent / 10. + 0.05);
+        probFunc[icent]->SetParameter(0, n_test);
+        probFuncB[icent]->SetParameter(0, sigma_Inelastic);
+        probFuncB[icent]->SetParameter(1, n_test);
+        for (int ipar = 0; ipar < 5; ipar++)
+        {
+            probFunc[icent]->SetParameter(ipar + 1, fitFunc->GetParameter(ipar));
+            probFuncB[icent]->SetParameter(ipar + 2, fitFunc->GetParameter(ipar));
         }
     }
 
@@ -393,14 +423,19 @@ int main(int argc, char **argv)
     fitFunc->Write();
     fitEtaFunc->Write();
     centFunc->Write();
-    for (int icent=0; icent<10; icent++)
+    for (int icent = 0; icent < 10; icent++)
     {
         probFunc[icent]->Write();
     }
-    for (int icent=0; icent<10; icent++)
+    for (int icent = 0; icent < 10; icent++)
     {
         probFuncB[icent]->Write();
     }
+    for (int icent = 0; icent < 10; icent++)
+    {
+        hBimp[icent]->Write();
+    }
+    hBimpTotal->Write();
 
     fo->Close();