From 59b46d22d6ab0607565497e7327465be03641d01 Mon Sep 17 00:00:00 2001 From: coppedis Date: Thu, 3 Jul 2014 12:59:14 +0200 Subject: [PATCH] Updated SNM Glauber fit --- .../GlauberSNM/AliCentralityGlauberFit.cxx | 581 ++++++++++++------ .../GlauberSNM/AliCentralityGlauberFit.h | 48 +- .../GlauberSNM/makeCentralityFitnozeri.C | 192 ++++++ 3 files changed, 624 insertions(+), 197 deletions(-) create mode 100644 PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C diff --git a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx index a95f2359560..3e003a4623e 100755 --- a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx +++ b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx @@ -34,7 +34,7 @@ #include #include #include -#include +#include #include #include #include @@ -57,9 +57,12 @@ AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) : fNbog(0), fBoglow(0), fBoghigh(0), - fNCP(0), - fCPlow(0), - fCPhigh(0), + fNgamma(0), + fgammalow(0), + fgammahigh(0), + fNsigmap(0), + fsigmaplow(0), + fsigmaphigh(0), fRebinFactor(0), fScalemin(0), fMultmin(0), @@ -73,13 +76,15 @@ AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) : fGlauberHist(0), fUseChi2(kTRUE), fNevents(100000), + fIsZN(kTRUE), + fNTrials(0), + fChi2Min(2000.), fInrootfilename(0), fInntuplename(0), fOutrootfilename(0), fOutntuplename(0), fAncfilename("ancestor_hists.root"), - fHistnames(), - fIsZN(kTRUE) + fHistnames() { // Standard constructor. TFile *f = 0; @@ -91,6 +96,7 @@ AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) : fGlauntuple->SetBranchAddress("B",&fB); fGlauntuple->SetBranchAddress("tAA",&fTaa); } + for(int i=0; i<4; i++) fParamnSlow[i] = 0.; } @@ -125,9 +131,12 @@ void AliCentralityGlauberFit::SetGlauberParam( Int_t Nbog, Double_t boglow, Double_t boghigh, - Int_t NCP, - Double_t CPlow, - Double_t CPhigh) + Int_t Ngamma, + Double_t gammalow, + Double_t gammahigh, + Int_t Nsigmap, + Double_t sigmaplow, + Double_t sigmaphigh) { // Set Glauber parameters. fNk=Nk; @@ -142,9 +151,23 @@ void AliCentralityGlauberFit::SetGlauberParam( fNbog=Nbog; fBoglow=boglow; fBoghigh=boghigh; - fNCP=NCP; - fCPlow=CPlow; - fCPhigh=CPhigh; + fNgamma=Ngamma; + fgammalow=gammalow; + fgammahigh=gammahigh; + fNsigmap=Nsigmap; + fsigmaplow=sigmaplow; + fsigmaphigh=sigmaphigh; +} + +//-------------------------------------------------------------------------------------------------- +void AliCentralityGlauberFit::SetNParam(Double_t a, Double_t b, Double_t c, Double_t d) +{ + // Setting parameters that are adjusted + fParamnSlow[0] = a; + fParamnSlow[1] = b; + fParamnSlow[2] = c; + fParamnSlow[3] = d; + for(int i=0; i<4; i++) printf(" INITIAL PARAMETER VALUES : paramnSlow[%d] = %f\n", i, fParamnSlow[i]); } //-------------------------------------------------------------------------------------------------- @@ -152,15 +175,14 @@ void AliCentralityGlauberFit::MakeFits() { // Make fits. - TH1F *hDATA; + TH1F *hDATA = 0x0; TH1F *thistGlau; TFile *inrootfile; TFile *outrootfile; - //FILE* fTxt = fopen ("parameters.txt","w"); // open inrootfile, outrootfile std::cout << "input file " << fInrootfilename << std::endl; - std::cout << "output file " << fOutrootfilename << std::endl; + std::cout << "output file " << fOutrootfilename << std::endl << std::endl; inrootfile = new TFile(fInrootfilename,"OPEN"); outrootfile = new TFile(fOutrootfilename,"RECREATE"); @@ -168,93 +190,132 @@ void AliCentralityGlauberFit::MakeFits() std::vector::const_iterator hni; for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) { hDATA = (TH1F*) (inrootfile->Get(*hni)); - std::cout << " -> Getting histogram " << *hni << std::endl; + std::cout << " -> Getting histogram " << *hni << std::endl << std::endl; if (!hDATA) { - TList *list = (TList*) (inrootfile->Get("CentralityStat")); - //TList *list = (TList*) (inrootfile->Get("VZEROEquaFactorStat")); - hDATA = (TH1F*) (list->FindObject(*hni)); + printf(" @@@@@@ No DATA histo in input -> EXITING!!!!!\n\n"); + return; + //TList *list = (TList*) (inrootfile->Get("CentralityStat")); + //hDATA = (TH1F*) (list->FindObject(*hni)); } //hDATA->Rebin(fRebinFactor); //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1)); - Double_t chi2min = 9999999.0; - Double_t CP_min=-1; + Double_t gamma_min=-1; + Double_t sigmap_min=-1; Double_t bog_min=-1; Double_t sigma_min=-1; Double_t alpha_min=-1; - Double_t k_min=-1; - Double_t alpha, k, sigma, bog, CP, chi2; - - for (Int_t nCP=0;nCP Chi2 %f \n", chi2); + if(chi2 < fChi2Min){ + fChi2Min = chi2; + alpha_min = alpha; + sigma_min = sigma; + bog_min = bog; + gamma_min = gamma; + sigmap_min = sigmap; + k_min = k; + } + } } } } } + } } - thistGlau->Reset(); - thistGlau = GlauberHisto(k_min,alpha_min,sigma_min,bog_min,CP_min, hDATA,kTRUE); - + if(fNsigmap>1 || fNgamma>1 || fNbog>1 || fNsigma>1 || fNalpha>1 || fNk>1 || fNTrials>1){ + thistGlau->Reset(); + printf(" **** calling GlauberHisto with acc = %1.3f R = %1.4f bog = %1.3f gamma = %1.3f sigmap = %1.3f\n", + alpha, sigma, bog, gamma, sigmap); + thistGlau = GlauberHisto(k_min,alpha_min,sigma_min,bog_min,gamma_min, sigmap_min, hDATA,kTRUE); + } + TH1F * hGLAU = 0x0; hGLAU = (TH1F *) thistGlau->Clone("hGLAU"); + hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU"))); - hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f_%.3f", - k_min,alpha_min,sigma_min,bog_min,CP_min))); + if(fIsZN) hGLAU->SetTitle(((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.2f-%.2f-%.3f-%.2f", + bog,sigmap,gamma,fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3]))); + else hGLAU->SetTitle(((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f", + alpha,bog,sigmap))); + Int_t firstBin=1; + if(fIsZN) firstBin = 8;//ZN Int_t lastBin=0; - if(fIsZN) lastBin = 100;//ZN + if(fIsZN) lastBin = 800;//ZN else lastBin = hDATA->GetNbinsX(); // ZP - Double_t mcintegral = hGLAU->Integral(2, lastBin,"width"); - Double_t dataintegral = hDATA->Integral(2, lastBin,"width"); - Double_t scale = (dataintegral/mcintegral); // - std::cout << "hGLAU -> Integral in bin range:" << "1-" << lastBin << " = " << - hGLAU->Integral(2, lastBin) << std::endl; - std::cout << "hDATA -> Integral in bin range:" << "1-" << lastBin << " = " << - hDATA->Integral(2, lastBin) << std::endl; + Double_t mcintegral = hGLAU->Integral(firstBin, lastBin,"width"); + Double_t dataintegral = hDATA->Integral(firstBin, lastBin,"width"); + Double_t scale = 1.; + if(mcintegral>0.) (scale = dataintegral/mcintegral); + printf("\n scale %f \n", scale); // - //printf(" binwidth: hGLAU %f hDATA %f\n", hGLAU->GetBinWidth(100), hDATA->GetBinWidth(100)); - printf(" scale %f \n", scale); + std::cout << "hGLAU -> Integral in bin range:" << firstBin << "-" << lastBin << " = " << + mcintegral << std::endl; + std::cout << "hDATA -> Integral in bin range:" << firstBin << "-" << lastBin << " = " << + dataintegral << std::endl; hGLAU->Scale(scale); // - std::cout << "hGLAU -> Integral (whole range): " << hGLAU->Integral(1, hGLAU->GetNbinsX()) << std::endl; - std::cout << "hDATA -> Integral (whole range): " << hDATA->Integral(1, hDATA->GetNbinsX()) << std::endl; + std::cout << "hGLAU -> Integral (whole range): " << hGLAU->Integral(1, hGLAU->GetNbinsX(),"width") << std::endl; + std::cout << "hDATA -> Integral (whole range): " << hDATA->Integral(1, hDATA->GetNbinsX(),"width") << std::endl; + + TCanvas *c2 = new TCanvas("c2","checkGLAU",100,0,700,700); + c2->cd(); + hGLAU->SetLineColor(kPink); + hGLAU->Draw(""); + hDATA->Draw("same"); SaveHisto(hDATA, hGLAU, outrootfile); - //fclose (fTxt); - printf("\n \t k = %1.2f alpha = %1.2f sigma = %1.2f bog = %1.2f CP = %1.2f \n\n", k,alpha, sigma, bog, CP); - std::cout << "chi2 min is " << chi2min << std::endl; + printf("\n \t k = %1.4f alpha = %1.3f sigma = %1.2f bog = %1.4f gamma = %1.4f sigmap = %.3f \n", k,alpha, sigma, bog, gamma, sigmap); + if(fIsZN) printf(" \t a = %1.2f b = %1.3f c = %1.2f d = %1.2f\n\n", fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3]); + std::cout << "chi2 min is " << fChi2Min << std::endl << std::endl; std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin), hGLAU->FindBin(fMultmax))/hGLAU->Integral() - << " of the total cross section" << std::endl; - fTempHist=hDATA; - fGlauberHist=hGLAU; + << " of the total cross section" << std::endl << std::endl; + fTempHist = hDATA; + fGlauberHist = hGLAU; } // close inrootfile, outrootfile @@ -264,84 +325,115 @@ void AliCentralityGlauberFit::MakeFits() //-------------------------------------------------------------------------------------------------- -TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog, Double_t CP, +TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog, Double_t gamma, Double_t sigmap, TH1F *hDATA, Bool_t save) { // Get Glauber histogram. - static TH1F *h1 = (TH1F*)hDATA->Clone("h1"); + static TH1F *h1 = (TH1F*) hDATA->Clone("h1"); h1->Reset(); h1->SetName(Form("fit_%.3f_%.3f",k,alpha)); - + TFile *outFile = NULL; TNtuple *ntuple = NULL; - if (save) { + if(save){ outFile = new TFile(fOutntuplename,"RECREATE"); - ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nbn:ngn:Etot"); + ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nblack:ngray:Etot:lcp:nslowp:nslown"); } - Int_t nents = 0; - if (fGlauntuple) - nents = fGlauntuple->GetEntries(); + Int_t nents=0, evtot=0, evn=0, evzeron=0, evzerop=0, evenzero=0; + if(fGlauntuple) nents = fGlauntuple->GetEntries(); for (Int_t i=0;iGetEntry(i % nents); - else { + if(fGlauntuple) fGlauntuple->GetEntry(i % nents); + else{ + printf(" NOT GETTING ENTRY FROM TNTUPLE for ev. %d -> setting Ncoll=1 \n\n",i); fNpart = 2; fNcoll = 1; } + //printf(" Getting entry %d from ntuple\n",i%nents); // Slow Nucleon Model from Chiara - Double_t ntot=0., n=0.; + Double_t ntot=0., n=0., lcp=0., np=0., nn=0.; Double_t nbn=0., ngn=0., nbp=0., ngp=0.; - MakeSlowNucleons2(fNcoll,alpha,k,bog,CP, nbn,ngn,nbp,ngp); - //MakeSlowNucleons2s(fNcoll,alpha,k,bog,CP, nbn,ngn,nbp,ngp); - - // acceptance + //MakeSlowNucleons(fNcoll,nbn,ngn,nbp,ngp); // pure Sikler + //MakeSlowNucleons2(fNcoll,bog,gamma, nbn,ngn,nbp,ngp,lcp); // smear in Ncoll + MakeSlowNucleons2s(fNcoll,bog,gamma,sigmap, nbn,ngn,nbp,ngp,lcp); // smear in Nslow + //MakeSlowNucleonsAM(fNcoll,bog,gamma, nbn,ngn,nbp,ngp); + + evtot++; + if(fNcoll==1) evn++; // - if(fIsZN) n = alpha*ngn+alpha*nbn; // ZNA - else n = alpha*ngp+alpha*nbp; // ZPA - //---------------------------------------- + if(fIsZN) n = alpha*(ngn+nbn); // ZNA + else n = alpha*(ngp+nbp); // ZPA + + np = nbp+ngp; + nn = nbn+ngn; + + if((alpha*(ngn+nbn)<0.5)) evzeron++; + if((alpha*(ngp+nbp)<0.5)) evzerop++; //------ experimental resolution -> Gaussian smearing ------------------------------------ - Double_t resexp=0; + Double_t resexp = 0.; //if (n>0) resexp = sigma*gRandom->Gaus(0.0,1.0)/TMath::Sqrt(n); //else resexp=0; //ntot = n*(1+resexp); - //if(n>0) resexp = sigma*TMath::Sqrt(n); - //else resexp=0; - ntot = (int) (gRandom->Gaus(n, resexp)); + if(n>=0){ + resexp = sigma*TMath::Sqrt(n); + ntot = (gRandom->Gaus(n, resexp)); + } + else{ + n=0.; + } + //if(ntot<1e-6) printf(" ntot %f -> n %f nbp %f ngp %f nbn %f ngn %f \n",ntot,n,nbp,ngp,nbn,ngn); //---------------------------------------- // non-lineary effect ------------------------------- - //ntot = k*ntot; - //ntot = ntot + k*ntot*ntot; + ntot = ntot + k*ntot*ntot; - //cout << ngn << " " << nbn << " " << ntot << endl; + Double_t nFact = 4.*82./208.; - Double_t nFact = 1.577; - Double_t Etot = nFact*ntot; + int ntote = (int) n; + Float_t Etot = nFact*ntote; + + //if(Etot<1e-6) printf(" Ncoll %1.2f -> n %f (%d) -> E %f nbp %f ngp %f nbn %f ngn %f \n", fNcoll,n,ntote,Etot,nbp,ngp,nbn,ngn); // + // USE ONLY IF ZDC TDC IS REQUESTED TO CONDITION THE SPECTRUM!!!!!!!!! + float resexpe=0.; + if(n>0.){ + resexpe = 1./TMath::Sqrt(n)*sigma*gRandom->Gaus(0.0, 1.0); + //resexpe = 1./TMath::Sqrt(n)*sigma; + // + Etot = Etot*(1+resexpe); + //Etot = gRandom->Gaus(Etot, resexpe); + // + // printf("resexpe = %f -> Etot = %f \n", resexpe, Etot); + } + // non-lineary effect ------------------------------- + Etot = Etot + k*Etot*Etot; - if (n>0) - resexp = 1./TMath::Sqrt(n)*sigma*gRandom->Gaus(0.0,1.0); - - Etot = Etot*(1+resexp); - //printf(" ntot %f Etot %f \n", ntot, Etot); - - if(ntot>0) { + if(Etot<0.5) evenzero++; + + if(Etot>0) { h1->Fill(Etot); - if (save) - ntuple->Fill(fNpart, fNcoll, fB, fTaa, ntot, nbn, ngn, Etot); + if(save) ntuple->Fill(fNpart, fNcoll, fB, fTaa, (float) ntot, nbn, ngn, Etot, lcp, np, nn); } } - if (save) { + if(save) { + //printf("\n ******** Fraction of events with Ncoll=1 in the ntuple %f\n",(float) (1.*evn/evtot)); + if(fIsZN){ + printf(" ******** Fraction of events with no neutrons %f\n",(float) (1.*evzeron/evtot)); + printf(" ******** Fraction of events with no neutron energy %f\n",(float) (1.*evenzero/evtot)); + printf(" DATA: fraction of events with no ZN signal = %f\n\n", 1-alpha); + } + else printf("\n ******** Fraction of events with no protons %f \n",(float) (1.*evzerop/evtot)); + ntuple->Write(); outFile->Close(); } + return h1; } @@ -351,8 +443,10 @@ Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau) // Note that for different values of neff the mc distribution is identical, just re-scaled below // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!! - Int_t lowchibin = hDATA->FindBin(fMultmin); - Int_t highchibin = hDATA->FindBin(fMultmax); + Int_t lowchibin = 8; + Int_t highchibin = 0; + if(fIsZN) highchibin = 800; + else highchibin = 25; Double_t mcintegral = thistGlau->Integral(1, thistGlau->GetNbinsX()); Double_t scale = (hDATA->Integral(1, hDATA->GetNbinsX())/mcintegral); @@ -406,61 +500,56 @@ void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TFile *outroot } //-------------------------------------------------------------------------------------------------- -void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Double_t k, Double_t &nbn, Double_t &ngn) +void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp) { // from AliGenSlowNucleonModelExp (Chiara Oppedisano) // Return the number of black and gray nucleons // // Number of collisions - Int_t fP = 82; - Int_t fN = 126; + Float_t fP = 82; + Float_t fN = 126; Float_t nu = (Float_t) ncoll; // // Mean number of gray nucleons - Float_t nGray = alpha * nu; + Float_t nGray = 2. * nu; Float_t nGrayNeutrons = nGray * fN / (fN + fP); Float_t nGrayProtons = nGray - nGrayNeutrons; // Mean number of black nucleons Float_t nBlack = 0.; - if(!fApplySaturation || (fApplySaturation && nGray=fnGraySaturation) nBlack = fnBlackSaturation; + if(nGray>=15.) nBlack = 28.; Float_t nBlackNeutrons = nBlack * 0.84; Float_t nBlackProtons = nBlack - nBlackNeutrons; -// Actual number (including fluctuations) from binomial distribution - Double_t p, ngp, nbp; - +// Actual number (including fluctuations) from binomial distribution // gray neutrons - p = nGrayNeutrons/fN; - ngn = gRandom->Binomial((Int_t) fN, p); + ngn = (double) gRandom->Binomial((Int_t) fN, nGrayNeutrons/fN); // gray protons - p = nGrayProtons/fP; - ngp = gRandom->Binomial((Int_t) fP, p); + ngp = (double) gRandom->Binomial((Int_t) fP, nGrayProtons/fP); // black neutrons - p = nBlackNeutrons/fN; - nbn = gRandom->Binomial((Int_t) fN, p); + nbn = (double) gRandom->Binomial((Int_t) fN, nBlackNeutrons/fN); // black protons - p = nBlackProtons/fP; - nbp = gRandom->Binomial((Int_t) fP, p); + nbp = (double) gRandom->Binomial((Int_t) fP, nBlackProtons/fP); + + // printf(" Ncoll %d ngp %f nbp %f ngn %f nbn %f \n",ncoll,ngp, nbp, ngn, nbn); } //-------------------------------------------------------------------------------------------------- - void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp) +void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Double_t gamma, + Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp, Double_t &lcp) { // from AliGenSlowNucleonModelExp (Chiara Oppedisano) // // Return the number of black and gray nucleons // -// Number of collisions - // based on E910 model ================================================================ + // PROTONS ---------------------------------------------------------------- Int_t fP = 82; Int_t fN = 126; @@ -468,10 +557,7 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub Float_t nu = (Float_t) ncoll; // nu = gRandom->Gaus(nu,0.5); - if(nu<1.) nu = 1.; - // - //float sdp = gRandom->Rndm(); - //if(nu==1. && sdp<=0.2) return; + //if(nu<1.) nu = 1.; // gray protons Float_t poverpd = 0.843; @@ -481,71 +567,198 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub Double_t p; p = nGrayp/fP; ngp = gRandom->Binomial((Int_t) fP, p); - //ngp = gRandom->Gaus(nGrayp, TMath::Sqrt(fP*p*(1-p))); - if(nGrayp<0.) ngp=0; + //if(nGrayp<0.) ngp=0; // black protons //Float_t blackovergray = 3./7.;// from spallation //Float_t blackovergray = 0.65; // from COSY Float_t blackovergray = bog; - // Float_t blackovergray = 0.65; // from COSY Float_t nBlackp = blackovergray*nGrayp; p = nBlackp/fP; nbp = gRandom->Binomial((Int_t) fP, p); - //nbp = gRandom->Gaus(nBlackp, TMath::Sqrt(fP*p*(1-p))); - if(nBlackp<0.) nbp=0; + // + //if(nBlackp<0.) nbp=0; + // NEUTRONS ---------------------------------------------------------------- + if(nu<3.){ nGrayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu; nBlackp = blackovergray*nGrayp; - } - + } // gray neutrons Float_t nGrayNeutrons = 0.; Float_t nBlackNeutrons = 0.; - Float_t cp = (nGrayp+nBlackp)/CP; + lcp = (nGrayp+nBlackp)/gamma; + + Float_t nSlow = 0.; + if(lcp>0.){ + //Float_t fParamnSlow[3] = {61., 470., 7.}; + nSlow = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+lcp)+fParamnSlow[3]*lcp; + float xconj = 3.; + float paramRetta = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+xconj); + if(lcp~0.50*) + //nGrayNeutrons = 0.47 * 2. * nu; // a la DPMJET (~0.39*) + nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (~0.11*) + } + // + if(nGrayNeutrons<0.) nGrayNeutrons=0.; + if(nBlackNeutrons<0.) nBlackNeutrons=0.; - // if(cp>0.){ - // //Float_t nSlow = 51.5+469.2/(-8.762-cp); - // Float_t nSlow = 60.0+469.2/(-8.762-cp); - // //if(cp<2.5) nSlow = 1+(9.9-1)/(2.5-0)*(cp-0); - // if(cp<3.) nSlow = 0.+(11.6-0.)/(3.-0.)*(cp-0.); - - // nGrayNeutrons = nSlow * 0.1; - // nBlackNeutrons = nSlow - nGrayNeutrons; - // } - if(cp>0.){ - Float_t paramnSlow[3] = {61., 470., 7.}; - Float_t nSlow = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-cp); - float paramRetta = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-3); - if(cp<3.) nSlow = 0.+(paramRetta-0.)/(3.-0.)*(cp-0.); + p = nGrayNeutrons/fN; + ngn = gRandom->Binomial((Int_t) fN, p); + + // black neutrons + p = nBlackNeutrons/fN; + //nbn = gRandom->Binomial((Int_t) fN, p); + nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p))); + if(nbn<0.) nbn=0.; + +} + +//-------------------------------------------------------------------------------------------------- + void AliCentralityGlauberFit::MakeSlowNucleons2s(Int_t ncoll, Double_t bog, Double_t gamma, Double_t sigmap, + Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp, Double_t &lcp) +{ +// from AliGenSlowNucleonModelExp (Chiara Oppedisano) +// +// Return the number of black and gray nucleons +// + + Int_t fP = 82; + Int_t fN = 126; + + Float_t nu = (Float_t) ncoll; + + // PROTONS ---------------------------------------------------------------- + // gray protons + Float_t poverpd = 0.843; + Float_t zAu2zPb = 82./79.; + Float_t grayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb; + //if(nu<3.) grayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu; + // + Float_t nGrayp = gRandom->Gaus(grayp, sigmap); +//changed +// Float_t nGrayp = gRandom->Gaus(grayp, (0.8*sigmap-sigmap)*grayp/90+sigmap); + if(nGrayp<0.) nGrayp=0; + + Double_t p=0.; + p = nGrayp/fP; + ngp = (double) gRandom->Binomial(fP, p); + + // black protons + //Float_t blackovergray = 3./7.;// =0.43 from spallation + //Float_t blackovergray = 0.65; // from COSY + Float_t blackovergray = bog; + // + Float_t nBlackp = blackovergray*nGrayp; + if(nBlackp<0.) nBlackp=0; + + p = nBlackp/fP; + nbp = (double) gRandom->Binomial(fP, p); + + // NEUTRONS ---------------------------------------------------------------- + // Relevant ONLY for n + // NB DOESN'T WORK FOR SMEARING IN Nslow!!!!!!!!!! + /*if(nu<=3.){ + grayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu; + // smearing + nGrayp = gRandom->Gaus(grayp, sigmap); + nBlackp = blackovergray*nGrayp; + }*/ - nGrayNeutrons = nSlow * 0.1; - nBlackNeutrons = nSlow - nGrayNeutrons; + // gray neutrons + lcp = (nGrayp+nBlackp)/gamma; + + Float_t nGrayNeutrons = 0.; + Float_t nBlackNeutrons = 0.; + Float_t nSlow = 0.; + if(lcp>0.){ + nSlow = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+lcp)+fParamnSlow[3]*lcp; + // +//changed +// float xconj = fParamnSlow[1]/fParamnSlow[0]-fParamnSlow[2]; + float xconj = 2.; + float yconj = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+xconj)+fParamnSlow[3]*xconj; + if(lcp ncoll = %1.0f -> ngrayn = %1.0f nblackn = %1.0f \n", nu, nGrayNeutrons, nBlackNeutrons); + // taking into account the prob4bility p to have 0 neutrons when we have 0 protons + // p = 1.-(0.82*0.96) ~ 0.22 (assuming that 100% of ev. without n are also without p) + //Float_t r = gRandom->Rndm(); + //if(r>0.1){ + // Sikler corrected for neutron yield... + nBlackNeutrons = 126./208. * 4. * nu; + //nGrayNeutrons = 126./208. * 2. * nu; // a la Sikler (~0.50*) + //nGrayNeutrons = 0.47 * 2. * nu; // a la DPMJET (~0.39*) + //nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (~0.11*) + // smearing! +//changed + nBlackNeutrons = gRandom->Gaus(nBlackNeutrons, sigmap); + //nGrayNeutrons = nBlackNeutrons/2.; // a la Sikler (~0.50*) + nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (~0.11*) + //printf(" Sikler -> Ncoll %f %f %f \n",nu,nGrayNeutrons, nBlackNeutrons); + //} } + // + if(nGrayNeutrons<0.) nGrayNeutrons=0.; + if(nBlackNeutrons<0.) nBlackNeutrons=0.; p = nGrayNeutrons/fN; //ngn = gRandom->Binomial((Int_t) fN, p); - ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN*p*(1-p))); +//changed + if(nGrayNeutrons>=10) ngn = gRandom->Gaus(nGrayNeutrons+0.5, TMath::Sqrt(fN*p*(1-p))); + else ngn = gRandom->PoissonD(nGrayNeutrons); + //else ngn = gRandom->Binomial((Int_t) fN, p); + if(ngn<0.) ngn=0.; // black neutrons p = nBlackNeutrons/fN; //nbn = gRandom->Binomial((Int_t) fN, p); - nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p))); +//changed + if(nBlackNeutrons>=10) nbn = gRandom->Gaus(nBlackNeutrons+0.5, TMath::Sqrt(fN*p*(1-p))); + else nbn = gRandom->PoissonD(nBlackNeutrons); + //else nbn = gRandom->Binomial((Int_t) fN, p); + if(nbn<0.) nbn=0.; + + //if(lcp<=0.) printf(" = %1.2f -> Ncoll = %1.2f -> %1.2f -> %1.0f gray n, %1.0f black n \n", + // (nGrayp+nBlackp), nu, nGrayNeutrons+nBlackNeutrons, ngn, nbn); + + +/*if(ncoll==1){ + printf(" Ncoll %d ", ncoll); + printf(" Ngp %1.0f Nbp %1.0f", ngp, nbp); + printf(" Ngn %1.0f Nbn %1.0f \n", ngn, nbn); +}*/ +/*if(ngn+nbn<0.005){ + printf(" Ncoll %d ", ncoll); + printf(" %f %f \t Ngp %1.0f Nbp %1.0f \t LCP %f\n", nGrayp, nBlackp, ngp, nbp,lcp); + printf(" %f Ngn %1.0f Nbn %1.0f \n", nSlow, ngn, nbn); +}*/ + } + //-------------------------------------------------------------------------------------------------- - void AliCentralityGlauberFit::MakeSlowNucleons2s(Int_t ncoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp) +/* void AliCentralityGlauberFit::MakeSlowNucleonsAM(Int_t ncoll, Double_t alpha, Double_t k, Double_t bog, + Double_t gamma, Double_t &nbn, Double_t &ngn, Double_t &nbp, Double_t &ngp) { // from AliGenSlowNucleonModelExp (Chiara Oppedisano) // @@ -566,36 +779,35 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub Float_t zAu2zPb = 82./79.; Float_t grayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb; Float_t nGrayp = gRandom->Gaus(grayp, sigmap); - if(nGrayp<0.) nGrayp=0.; + //if(nGrayp<0.) nGrayp=0.; Double_t p=0.; p = nGrayp/fP; ngp = gRandom->Binomial((Int_t) fP, p); //ngp = gRandom->Gaus(nGrayp, TMath::Sqrt(fP*p*(1-p))); - if(nGrayp<0.) ngp=0; // black protons //Float_t blackovergray = 3./7.;// from spallation - Float_t blackovergray = 0.65; // from COSY - //Float_t blackovergray = bog; - Float_t blackp = blackovergray*nGrayp; - Float_t nBlackp = gRandom->Gaus(blackp, sigmap); - if(nBlackp<0.) nBlackp=0.; + //Float_t blackovergray = 0.65; // from COSY + Float_t blackovergray = bog; + // + //Float_t blackp = blackovergray*nGrayp; + //Float_t nBlackp = gRandom->Gaus(blackp, sigmap); + Float_t nBlackp = blackovergray*nGrayp; + //if(nBlackp<0.) nBlackp=0.; p = nBlackp/fP; nbp = gRandom->Binomial((Int_t) fP, p); - //nbp = gRandom->Gaus(nBlackp, TMath::Sqrt(fP*p*(1-p))); - if(nBlackp<0.) nbp=0; // gray neutrons Float_t nGrayNeutrons = 0.; Float_t nBlackNeutrons = 0.; - Float_t cp = (nGrayp+nBlackp)/CP; + Float_t cp = (nGrayp+nBlackp)/gamma; if(cp>0.){ - Float_t paramnSlow[3] = {61., 470., 7.}; - Float_t nSlow = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-cp); - //float paramRetta = paramnSlow[0]+paramnSlow[1]/(-paramnSlow[2]-3); + //Float_t fParamnSlow[3] = {61., 470., 7.}; + Float_t nSlow = fParamnSlow[0]+fParamnSlow[1]/(-fParamnSlow[2]-cp); + //float paramRetta = fParamnSlow[0]+fParamnSlow[1]/(-fParamnSlow[2]-3); //if(cp<3.) nSlow = 0.+(paramRetta-0.)/(3.-0.)*(cp-0.); nGrayNeutrons = nSlow * 0.1; @@ -603,8 +815,8 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub } else{ // Sikler "pasturato" - nGrayNeutrons = 0.47 * 2.2 * nu; // fAlphaGray=2.3 - nBlackNeutrons = 0.88 * 3.6 * nu; // fAlphaBlack=3.6 + nGrayNeutrons = 0.47 * 2.2 * nu; + nBlackNeutrons = 0.88 * 3.6 * nu; //printf("nslowp=0 -> ncoll = %1.0f -> ngrayn = %1.0f nblackn = %1.0f \n", nu, nGrayNeutrons, nBlackNeutrons); } // @@ -614,14 +826,21 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub p = nGrayNeutrons/fN; //ngn = gRandom->Binomial((Int_t) fN, p); ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN*p*(1-p))); - if(p<0.) ngn=0; // black neutrons p = nBlackNeutrons/fN; //nbn = gRandom->Binomial((Int_t) fN, p); nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p))); - if(p<0.) nbn=0; -} + + if(nu>12){ + //Float_t spread[20] = {0.98,1.99,2.74,3.53,4.22,4.77,5.46,5.90,6.36, + // 6.63,6.91,7.42,7.74,7.96,7.03,6.76,9.49,11.91,12.16,5.66}; + //int index = (int) (nu-13); + Double_t n = gRandom->Gaus(40.34, 10.); + ngn = 0.1*n; + nbn = n-ngn; + } +}*/ Double_t AliCentralityGlauberFit::ConvertToEnergy(Double_t T) { @@ -629,7 +848,7 @@ Double_t AliCentralityGlauberFit::ConvertToEnergy(Double_t T) const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass(); Double_t fPmax =10; - Double_t m, p, pmax, ekin, f; + Double_t m, p, pmax, f; m = kMassNeutron; pmax = TMath::Sqrt(2*T*(T+TMath::Sqrt(T*T+m*m))); diff --git a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h index 31fafe6f8fe..dcda9f95314 100755 --- a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h +++ b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h @@ -28,7 +28,10 @@ class AliCentralityGlauberFit : public TObject { void AddHisto(TString name) { fHistnames.push_back(name); } TH1F *GetGlauberHist() { return fGlauberHist; } void MakeFits(); - void SetGlauberParam(Int_t Nk, Double_t klow, Double_t khigh, Int_t Nalpha, Double_t alphalow, Double_t alphahigh, Int_t Nsigma, Double_t sigmalow, Double_t sigmahigh, Int_t Nbog, Double_t boglow, Double_t boghigh, Int_t NCP, Double_t CPlow, Double_t CPhigh); + void SetGlauberParam(Int_t Nk, Double_t klow, Double_t khigh, Int_t Nalpha, Double_t alphalow, Double_t alphahigh, + Int_t Nsigma, Double_t sigmalow, Double_t sigmahigh, Int_t Nbog, Double_t boglow, Double_t boghigh, + Int_t Ngamma, Double_t gammalow, Double_t gammahigh, Int_t Nsigmap, Double_t sigmaplow, Double_t sigmaphigh); + void SetNParam(double a, double b, double c, double d); void SetInputFile(TString filename) { fInrootfilename = filename; } void SetInputNtuple(TString ntuplename) { fInntuplename = ntuplename; } void SetNevents(Int_t f) { fNevents=f; } @@ -42,6 +45,8 @@ class AliCentralityGlauberFit : public TObject { void SetSaturationParams(Int_t ngray=15, Int_t nblack=28) {fnGraySaturation=ngray; fnBlackSaturation=nblack;} void SetIsZN() {fIsZN=kTRUE;} void SetIsZP() {fIsZN=kFALSE;} + void SetNTrials(int ntrials) {fNTrials=ntrials;} + void SetChi2Min(double chi2min) {fChi2Min=chi2min;} private: Int_t fNk; // k points @@ -53,12 +58,15 @@ class AliCentralityGlauberFit : public TObject { Int_t fNsigma; // sigma points Double_t fSigmalow; // sigma low Double_t fSigmahigh; // sigma high - Int_t fNbog; // bog points - Double_t fBoglow; // bog low - Double_t fBoghigh; // bog high - Int_t fNCP; // CP points - Double_t fCPlow; // CP low - Double_t fCPhigh; // CP high + Int_t fNbog; // bog points + Double_t fBoglow; // bog low + Double_t fBoghigh; // bog high + Int_t fNgamma; // gamma points + Double_t fgammalow; // gamma low + Double_t fgammahigh; // gamma high + Int_t fNsigmap; // sigmap points + Double_t fsigmaplow; // sigmap low + Double_t fsigmaphigh; // sigmap high Int_t fRebinFactor; // rebin factor Double_t fScalemin; // mult min where to scale Double_t fMultmin; // mult min @@ -75,7 +83,11 @@ class AliCentralityGlauberFit : public TObject { Bool_t fApplySaturation;// Int_t fnGraySaturation; // Int_t fnBlackSaturation;// - Bool_t fIsZN; + Bool_t fIsZN; + Float_t fSDprobability[40]; + Double_t fParamnSlow[4]; + Int_t fNTrials; + Double_t fChi2Min; TString fInrootfilename; // input root file TString fInntuplename; // input Gauber ntuple @@ -84,19 +96,23 @@ class AliCentralityGlauberFit : public TObject { TString fAncfilename; // ancestor file name std::vector fHistnames; // histogram names - Double_t CalculateChi2(TH1F *hDATA, TH1F *thistGlau); - TH1F *GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog, Double_t CP, TH1F *hDATA, Bool_t save=kFALSE); - void SaveHisto(TH1F *hist1,TH1F *hist2, TFile *outrootfile); - void MakeSlowNucleons(Int_t fNcoll, Double_t alpha, Double_t k, Double_t &nbn, Double_t &ngn); - void MakeSlowNucleons2(Int_t fNcoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp); - void MakeSlowNucleons2s(Int_t fNcoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp); - //void MakeSlowNucleons2(Int_t fNcoll, Double_t alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn); + Double_t CalculateChi2(TH1F *hDATA, TH1F *thistGlau); + TH1F *GlauberHisto(Double_t k, Double_t alpha, Double_t sigma, Double_t bog, Double_t gamma, Double_t sigmap, + TH1F *hDATA, Bool_t save=kFALSE); + void SaveHisto(TH1F *hist1,TH1F *hist2, TFile *outrootfile); + void MakeSlowNucleons(Int_t fNcoll, Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp); + void MakeSlowNucleons2(Int_t fNcoll, Double_t bog, Double_t gamma, + Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp, Double_t &lcp); + void MakeSlowNucleons2s(Int_t fNcoll, Double_t bog, Double_t gamma, Double_t sigmap, + Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp, Double_t &lcp); + //void MakeSlowNucleonsAM(Int_t fNcoll, Double_t bog, Double_t gamma, + // Double_t &nbn, Double_t &ngn,Double_t &nbp, Double_t &ngp); Double_t ConvertToEnergy(Double_t T); Double_t Maxwell(Double_t m, Double_t p, Double_t T); AliCentralityGlauberFit(const AliCentralityGlauberFit&); AliCentralityGlauberFit &operator=(const AliCentralityGlauberFit&); - ClassDef(AliCentralityGlauberFit, 2) + ClassDef(AliCentralityGlauberFit, 4) }; #endif diff --git a/PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C b/PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C new file mode 100644 index 00000000000..b298082e6ad --- /dev/null +++ b/PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C @@ -0,0 +1,192 @@ +#if !defined(__CINT__) || defined(__MAKECINT__) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliCentralityGlauberFit.h" + +#endif + +void makeCentralityFitnozeri(const int nRun=195483, const char *system = "ZNA", + int ntrials=1, int Rebin=1, int Nevt=1.e7) +{ + //load libraries + gSystem->SetBuildDir("."); + gSystem->Load("libCore.so"); + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libSTEERBase.so"); + gROOT->ProcessLine(".include $ALICE_ROOT/include"); + gROOT->LoadMacro("AliCentralityGlauberFit.cxx++g"); + + const char *finnameGlau ="GlauberMC_pPb_ntuple_sigma70_mind4_r662_a546_Rpro6.root"; + char histname[8]; + double chi2min=0.; + + if((strncmp(system,"ZNA",3))==0){ + printf("\n Glauber fit on ZNA spectrum\n\n"); + sprintf(histname,"hZNA"); + chi2min = 200.; + } + else if((strncmp (system,"ZPA",3)) == 0){ + printf("\n Glauber fit on ZPA spectrum\n\n"); + sprintf(histname,"hZPA"); + chi2min=5000.; + } + + TString finname = Form("centrHistos%d-dNdeta.root",nRun); + printf(" Opening file %s\n",finname.Data()); + // + TString foutname = Form("%s_fit_%d.root",system,nRun); + TString foutnameGlau = Form("%s_ntuple_%d.root",system,nRun); + + + AliCentralityGlauberFit *mPM = new AliCentralityGlauberFit(finnameGlau); + mPM->SetInputFile(finname); + mPM->SetInputNtuple(finnameGlau); + mPM->SetOutputFile(foutname); + mPM->SetOutputNtuple(foutnameGlau); + mPM->AddHisto(histname); + + mPM->SetRebin(Rebin); + mPM->SetNevents(Nevt); + mPM->UseChi2(kTRUE); // If TRUE minimize Chi2 + mPM->SetNTrials(ntrials); + mPM->SetChi2Min(chi2min); + // + // COSY ORIGINAL!!!! + //mPM->SetNParam(51.5, 469.2, 8.762); + // COSY-like +// mPM->SetNParam(51.5, 620., 10.); +// change 2 + mPM->SetNParam(29., 132., 4.46, 0.29); + //mPM->SetNParam(29., 133., 4.48, 0.42); + //else mPM->SetNParam(51.5, 580., 9.5); + // + // ALICE ex-novo + //mPM->SetNParam(78., 700., 7.); + //mPM->SetNParam(61., 470., 7.); + + if(strncmp(system,"ZNA",3) == 0) { + printf(" Setting parameters for ZNA Glauber fit\n\n"); + mPM->SetIsZN(); + mPM->SetRangeToFit(0.5, 90.5); + mPM->SetRangeToScale(0); + if(nRun==195483) mPM->SetGlauberParam(1,0.0,0.0, 1,0.957,1., 1,0.29,0.30, 1,0.65,0.65, 1,0.575,0.585, 1, 0.26,0.3); + else mPM->SetGlauberParam(1,0.0,0.0, 1,0.957,1., 1,0.25,0.30, 1,0.65,0.65, 1,0.56,0.585, 1, 0.26,0.3); + } + else if (strncmp(system,"ZPA",3) == 0) { + mPM->SetIsZP(); + mPM->SetRangeToFit(1., 30.); + mPM->SetRangeToScale(0); + if(nRun==195483) mPM->SetGlauberParam(1,0.0016,0.0025, 1,0.62,0.8, 1,0.5,0.65, 1,0.65,0.65, 1,0.56,0.585, 1, 0.26,0.4); + else mPM->SetGlauberParam(1,0.0000,0.0000, 1,0.60,0.8, 1,0.5,0.65, 1,0.65,0.65, 1,0.56,0.585, 1, 0.26,0.4); +// mPM->SetGlauberParam(1,0.,1., 1,0.8,0.9, 1,0.6,0.65, 1,0.65,0.65, 1,0.56,0.585, 1, 0.28,0.4); + } + mPM->MakeFits(); + + char hnam[10], hnamg[20]; + double xt=0; + if(strncmp(system,"ZNA",3) == 0){ + sprintf(hnam,"hZNA"); + sprintf(hnamg,"hZNA_GLAU"); + xt=20.; + } + else if(strncmp(system,"ZPA",3) == 0) { + sprintf(hnam,"hZPA"); + sprintf(hnamg,"hZPA_GLAU"); + xt=5.; + } + + TFile *f = TFile::Open(foutname); + TH1 * hd = dynamic_cast (f->Get((hnam))); + TH1 * hg = dynamic_cast (f->Get((hnamg))); + hg->SetLineColor(kPink-2); + hg->SetLineWidth(2); + hd->SetMarkerStyle(20); + //hd->SetMarkerSize(1.2); + hd->SetMarkerColor(kBlue+3); + hd->SetLineWidth(2); + //hd->SetMinimum(10.); + + TCanvas *g = new TCanvas("g","g",0,0,700,700); + g->cd(); + gPad->SetLogy(1); + hg->Draw("E"); + hd->Draw("PEsame"); + hd->SetXTitle("E_{ZNA} (TeV)"); + hg->Draw("Esame"); + + printf("\n Entries with zero energy %1.0f -> %f of the total (%d)\n\n", + hg->GetBinContent(hg->FindBin(0.)), hg->GetBinContent(hg->FindBin(0.))/hg->Integral(), + hg->Integral()); + + /*double zncut[7] = {0., 16.3423, 38.2513, 52.7698, 63.7601, 70.2207, 75.0126}; + TH1F *hist[7]; + for(int ih=0; ih<7; ih++){ + char hnam[24]; + sprintf(hnam,"hZNA%d",ih); + hist[ih] = new TH1F(hnam,hnam,hd->GetNbinsX(),0.,142.5); + } + int index=0; + for(int ic=0; ic<7; ic++){ + for(int ib=1; ibGetNbinsX(); ib++){ + if(hd->GetBinLowEdge(ib)>=zncut[ic] && hd->GetBinLowEdge(ib+1)SetBinContent(ib, hd->GetBinContent(ib)); + if(hd->GetBinError(ib)>0) hist[index]->SetBinError(ib, 1./(hd->GetBinError(ib)*hd->GetBinError(ib))); + } + } + } + for(int ih=0; ih<7; ih++){ + if(ih%2==0) hist[ih]->SetFillColor(kAzure+6); + hist[ih]->Draw("hist SAME"); + } + hg->Draw("Esame");*/ + + TLatex text0; + text0.SetTextSize(0.04); + text0.SetTextColor(kBlue+3); + char ch[60]; + if(strncmp(system,"ZNA",3) == 0) sprintf(ch," DATA = %1.2f TeV ",hd->GetMean()); + else sprintf(ch," DATA = %1.2f TeV ",hd->GetMean()); + text0.DrawLatex(xt,20.,ch); + char chd[60]; + if(strncmp(system,"ZNA",3) == 0) sprintf(chd," Glauber = %1.2f TeV",hg->GetMean()); + else sprintf(chd," Glauber = %1.2f TeV",hg->GetMean()); + text0.SetTextColor(kPink-2); + text0.DrawLatex(xt,10.,chd); + + char ct[60]; + sprintf(ct,"RUN %d",nRun); + text0.SetTextColor(kTeal+2); + text0.DrawLatex(xt+80.,1.e4,ct); + + char psn[30]; + if(strncmp(system,"ZNA",3) == 0) sprintf(psn,"ZNA_fit%d.gif",nRun); + else sprintf(psn,"ZPA_fit%d.gif",nRun); + g->Print(psn); + +} + -- 2.43.0