#include "CentralityEstimator.h" #include "TTree.h" #include "TCanvas.h" #include "TF1.h" #include "TFile.h" #include "TGraph.h" #include "TLegend.h" #include "TMath.h" #include "TRandom.h" ClassImp(CentralityEstimator) // ----- Default constructor ------------------------------------------- CentralityEstimator::CentralityEstimator() : TNamed(), fMode("Default") { } void CentralityEstimator::LoadInputTree(Int_t n) { fSim = new TFile (fInputfilename.Data()); if (!fSim) { std::cout << "LoadInputTree: *** No File: "<< fInputfilename <<" - Error - " << std::endl; exit(EXIT_FAILURE); } fSimTree = (TTree*) fSim->Get(fTreeName.Data()); if (!fSimTree) { std::cout << "LoadInputTree: *** No Tree:"<< fTreeName <<" - Error - " << std::endl; exit(EXIT_FAILURE); } fSimTree->SetBranchAddress("B", &fB); fSimTree->SetBranchAddress("Npart", &fNpart); fSimTree->SetBranchAddress("Ncoll", &fNcoll); if (n < 0){ nEntries = fSimTree->GetEntries(); std::cout << "SetInputTree: *** Warning - Number of Entries < 0 or not set. " "total event number "<< nEntries <<" will be used" << std::endl; } } Double_t CentralityEstimator::GetEstimator(Double_t f, TString flag) { Double_t Na=0; //if (f == 1) Na = fNpart; if (flag == "") Na = fNpart; else if (flag == "NpartThreshold") Na = (f + fNpart/fNpartMax) > 1. ? fNpart : fNpart * (1 - fNpart/fNpartMax ); else if (flag == "Default") Na = fNpart; else if (flag == "Scaling") Na = fNpart - f * (fNpartMax/fNpart); else if (flag == "NpartOnly") Na = f*fNpart + (1-f)*fNcoll; else if (flag == "Npart") Na = TMath::Power(fNpart, f); else if (flag == "Ncoll") Na = TMath::Power(fNcoll, f); else if (flag == "PSD") Na = f - fNpart; // // no lininiarity // Double_t p0 = 0.740718;// +/- 0.0264785 // Double_t p1 = 0.000187403;// +/- 0.000252608 // Double_t p2 = -8.92373e-07;// +/- 5.56444e-07 // non-liniarity - due reco efficiency // Double_t p0 = 0.784036;// +/- 0.0264785 // Double_t p1 = -1.16025e-06;// +/- 0.000252608 Double_t p0 = 1; Float_t nHits = 0; for (Int_t j=0; jRndm() < p0 + f *fNpart*fNpart ) nHits += (int)hNBD->GetRandom(); // if(gRandom->Rndm() < p0 + p1 *fNpart + p2 *fNpart*fNpart ) nHits += (int)hNBD->GetRandom(); // nHits += (int)hNBD->GetRandom(); // if(gRandom->Rndm() + f < fNpart/fNpartMax) nHits += (int)hNBD->GetRandom(); // nHits += GetRandomNBD(k, mu); } // std::cout << "fNpartMax " << fNpartMax << " nHits "<< nHits << " Npart X " << fNpart << " Na " << Na << " fMode "<< fMode << std::endl; return nHits; } void CentralityEstimator::AddEstimator(TString name, Double_t f, Double_t mu, Double_t k, TString mode) { cout << "Adding Estimator: " << name <<" f: " << f << " mu: " << mu << " k: " << k << " mode:" << mode << endl; estimatorName.push_back(name); estimatorMode.push_back(mode); estimatorF.push_back (f); estimatorMu.push_back (mu); estimatorK.push_back (k); } void CentralityEstimator::Run() { fNewfile = new TFile(fOutputfilename ,"recreate"); fNewtree = fSimTree->CloneTree(); for(UInt_t iArray=0;iArrayBranch( estimatorName[iArray].Data() ,&estimator ,Form("%s/F",estimatorName[iArray].Data()) ); SetNBDhist(estimatorMu[iArray], estimatorK[iArray]); for (Int_t i=0; iGetEntry(i); estimator = GetEstimator(estimatorF[iArray],estimatorMode[iArray]); bpt->Fill(); if (!(i%100)) cout << "Filling Estimator: " << estimatorName[iArray].Data() <<" Event # " << i << "\r" << flush; } cout << endl; } /* AddEstimator("TOF" , -1.64069e-06 , 0.202237 , 6.36083 ); //f = -1.64069e-06 mu = 0.202237 k = 6.36083 Chi2Min = 50.1303 AddEstimator("RPC" , -1.64069e-06 , 0.503523 , 29.063 ); // f = -1.64069e-06 mu = 0.503523 k = 29.063 Chi2Min = 39.0985 AddEstimator("TOFRPCnew" , 1. , 0.742696 , 20.5058 ); // f = 1 mu = 0.742696 k = 20.5058 Chi2Min = 1.03845 // AddEstimator("SelectedParticleCand" , 0.67 , 1.52855 , 908.567 ); // // f = 0.805 mu = 0.686607 k = 71.2954 Chi2Min = 91.9057 // // f = 0.67 mu = 1.52855 k = 908.567 sigma = 2.13642 chi2 = 30.6328 BestChi2= 30.6328 // // AddEstimator("SelectedParticleCand1" , 1. , 0.214634 , 9.33643 ); // //f = 1 mu = 0.214634 k = 9.33643 Chi2Min = 74.3104 // AddEstimator("PrimaryParticleCand" , -1.3825e-07 , 0.185581 , 7.60447 ); //f = -1.3825e-07 mu = 0.185581 k = 7.60447 Chi2Min = 89.24 //AddEstimator("PrimaryParticleCand" ,-3.31154e-07, 0.196316, 47.9113 ); //-3.31154e-07, 0.196316, 47.9113 // AddEstimator("PrimaryParticleCand" , -1.0975e-07 , 0.235526 , 20.3434 ); // f = -1.0975e-07 mu = 0.235526 k = 20.3434 Chi2Min = 97.6704 // AddEstimator("MdcWires" , 1, 0.377, 0.791 ); //no fit // AddEstimator("FWSumChargeSpec" , 1, 0.217, 1.151 ); AddEstimator("SelectedParticleCand" , -1.3825e-07 ,0.215636 , 5.18754 ); //f = -1.3825e-07 mu = 0.215636 k = 5.18754 Chi2Min = 45.4244 */ fNewtree->Print(); fNewtree->Write("", TObject::kOverwrite); fNewfile->Close(); } void CentralityEstimator::SetNBDhist(Double_t mu, Double_t k) { // Interface for TH1F. Int_t nBins = 300; hNBD = new TH1F("hNBD","",nBins, 0, nBins); hNBD->SetName(Form("nbd_%f_%f",mu,k)); hNBD->SetDirectory(0); for (Int_t i=0; i1e-20) hNBD->SetBinContent(i+1, val); // std::cout << "val " << val << std::endl; } } Double_t CentralityEstimator::NBD(Double_t n, Double_t mu, Double_t k) const { // Compute NBD. Double_t F; Double_t f; if (n+k > 100.0) { // log method for handling large numbers F = TMath::LnGamma(n + k)- TMath::LnGamma(n + 1.)- TMath::LnGamma(k); f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k); F = F+f; F = TMath::Exp(F); } else { F = TMath::Gamma(n + k) / ( TMath::Gamma(n + 1.) * TMath::Gamma(k) ); f = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k); f = TMath::Exp(f); F *= f; } return F; } Int_t CentralityEstimator::GetRandomNBD(Double_t k, Double_t nbar) const { // negative binomial distribution generator, S. Voloshin, 09-May-2007 Double_t sum=0.; Int_t i=0; Double_t ran=gRandom->Rndm(); Double_t trm=1./pow(1.+nbar/k,k); if (trm==0.) { cout<<"NBD overflow"<<" nbar="<