#include <TGraphErrors.h>
#include <vector>
#include <TMinuit.h>
-#include <TStopwatch.h>
+#include <TCanvas.h>
#include <TRandom.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
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),
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;
fGlauntuple->SetBranchAddress("B",&fB);
fGlauntuple->SetBranchAddress("tAA",&fTaa);
}
+ for(int i=0; i<4; i++) fParamnSlow[i] = 0.;
}
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;
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]);
}
//--------------------------------------------------------------------------------------------------
{
// 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");
std::vector<TString>::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<fNCP;nCP++) {
- CP = fCPlow + ((Double_t) nCP ) * (fCPhigh - fCPlow ) / fNCP;
+ Double_t k_min=-1, chi2=0.;
+ Double_t alpha=0., k=0., sigma=0., bog=0., gamma=0., sigmap=0.;
+
+ for(Int_t nsigmap=0; nsigmap<fNsigmap; nsigmap++) {
+ sigmap = fsigmaplow + ((Double_t) nsigmap ) * (fsigmaphigh - fsigmaplow ) / fNsigmap;
+
+ for(Int_t ngamma=0; ngamma<fNgamma; ngamma++) {
+ gamma = fgammalow + ((Double_t) ngamma ) * (fgammahigh - fgammalow ) / fNgamma;
- for (Int_t nbog=0;nbog<fNbog;nbog++) {
+ for(Int_t nbog=0; nbog<fNbog; nbog++) {
bog = fBoglow + ((Double_t) nbog ) * (fBoghigh - fBoglow ) / fNbog;
- for (Int_t nsigma=0;nsigma<fNsigma;nsigma++) {
+ for(Int_t nsigma=0; nsigma<fNsigma; nsigma++) {
sigma = fSigmalow + ((Double_t) nsigma ) * (fSigmahigh - fSigmalow ) / fNsigma;
- for (Int_t nalpha=0;nalpha<fNalpha;nalpha++) {
+ for(Int_t nalpha=0; nalpha<fNalpha; nalpha++) {
alpha = fAlphalow + ((Double_t) nalpha ) * (fAlphahigh - fAlphalow ) / fNalpha;
- for (Int_t nk=0; nk<fNk; nk++) {
+ for(Int_t nk=0; nk<fNk; nk++) {
k = fKlow + ((Double_t) nk ) * (fKhigh - fKlow ) / fNk;
-
- thistGlau = GlauberHisto(k,alpha,sigma,bog,CP,hDATA,kFALSE);
- chi2 = CalculateChi2(hDATA,thistGlau);
- //cout << "chi2 " << chi2<< endl;
- if ( chi2 < chi2min ) {
- chi2min = chi2;
- alpha_min=alpha;
- sigma_min=sigma;
- bog_min=bog;
- CP_min=CP;
- k_min=k;
- }
+
+
+ for(int i=0; i<fNTrials; i++){
+ printf(" **** calling GlauberHisto with acc = %1.3f R = %1.4f k = %1.4f bog = %1.3f gamma = %1.3f sigmap = %1.3f\n",
+ alpha, sigma, k, bog, gamma, sigmap);
+ if(fIsZN){
+ // Updating parameters for slow n
+ fParamnSlow[0] = fParamnSlow[0]; // + 0.5*i;
+ fParamnSlow[1] = fParamnSlow[1]; // + 10.*i;
+ fParamnSlow[2] = fParamnSlow[2] + 0.05*i;
+ fParamnSlow[3] = fParamnSlow[3];// + 0.01*i;
+ //
+ printf(" \t\t trial #%d n slow param.: %1.2f %1.2f %1.2f %1.2f\n",
+ i,fParamnSlow[0], fParamnSlow[1], fParamnSlow[2], fParamnSlow[3]);
+ }
+
+ //thistGlau = GlauberHisto(k,alpha,sigma,bog,gamma,sigmap,hDATA,kFALSE);
+ thistGlau = GlauberHisto(k,alpha,sigma,bog,gamma,sigmap,hDATA,kTRUE);
+ //
+ chi2 = CalculateChi2(hDATA, thistGlau);
+ printf(" -> 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
//--------------------------------------------------------------------------------------------------
-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;i<fNevents;++i) {
- if (fGlauntuple)
- fGlauntuple->GetEntry(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;
}
// 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);
}
//--------------------------------------------------------------------------------------------------
-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 = k * nu;
- else if(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;
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;
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<xconj) nSlow = 0.+(paramRetta-0.)/(xconj-0.)*(lcp-0.);
+ //
+ nSlow = 1.5*nSlow;
+ //
+ nBlackNeutrons = nSlow * 0.9;;
+ nGrayNeutrons = nSlow - nBlackNeutrons;
+ }
+ else{
+ // slightly adjusted Sikler corrected for neutron yield...
+ nBlackNeutrons = 126./208. * 3.6 * nu;
+ //nGrayNeutrons = 126./208. * 2. * nu; // a la Sikler (<Ngn>~0.50*<Nbn>)
+ //nGrayNeutrons = 0.47 * 2. * nu; // a la DPMJET (<Ngn>~0.39*<Nbn>)
+ nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (<Ngn>~0.11*<Nbn>)
+ }
+ //
+ 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<xconj) nSlow = 0.+(yconj-0.)*(lcp-0.)/(xconj-0.);
+ //
+ // Factor to scale COSY to obtain ALICE n mltiplicities
+ nSlow = 1.68*nSlow;
+ //
+ nBlackNeutrons = nSlow * 0.9;;
+ nGrayNeutrons = nSlow - nBlackNeutrons;
}
else{
- // Sikler "pasturato"
- nGrayNeutrons = 0.47 * 2.2 * nu; // fAlphaGray=2.3
- nBlackNeutrons = 0.88 * 3.6 * nu; // fAlphaBlack=3.6
- if(nGrayNeutrons<0.) nGrayNeutrons=0.;
- if(nBlackNeutrons<0.) nBlackNeutrons=0.;
- //printf("nslowp=0 -> 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 (<Ngn>~0.50*<Nbn>)
+ //nGrayNeutrons = 0.47 * 2. * nu; // a la DPMJET (<Ngn>~0.39*<Nbn>)
+ //nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (<Ngn>~0.11*<Nbn>)
+ // smearing!
+//changed
+ nBlackNeutrons = gRandom->Gaus(nBlackNeutrons, sigmap);
+ //nGrayNeutrons = nBlackNeutrons/2.; // a la Sikler (<Ngn>~0.50*<Nbn>)
+ nGrayNeutrons = nBlackNeutrons/9.; // a la spallation (<Ngn>~0.11*<Nbn>)
+ //printf(" Sikler -> Ncoll %f <Ngrayn> %f <Nblackn> %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("<Nslowp> = %1.2f -> Ncoll = %1.2f -> <Nslown> %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(" <Ngp> %f <Nbp> %f \t Ngp %1.0f Nbp %1.0f \t LCP %f\n", nGrayp, nBlackp, ngp, nbp,lcp);
+ printf(" <Nslow,n> %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)
//
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;
}
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);
}
//
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)
{
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)));
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include <cstdlib>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TROOT.h>
+#include <Riostream.h>
+#include <TSystem.h>
+#include <TClassTable.h>
+#include <TStyle.h>
+#include <TMath.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TH1F.h>
+#include <TH1D.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TLine.h>
+#include <TNtuple.h>
+#include <TLatex.h>
+#include <TGraphErrors.h>
+#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<TH1*> (f->Get((hnam)));
+ TH1 * hg = dynamic_cast<TH1*> (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; ib<hd->GetNbinsX(); ib++){
+ if(hd->GetBinLowEdge(ib)>=zncut[ic] && hd->GetBinLowEdge(ib+1)<zncut[ic+1]){
+ index = ic;
+ hist[index]->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,"<E_{ZNA}> DATA = %1.2f TeV ",hd->GetMean());
+ else sprintf(ch,"<E_{ZPA}> DATA = %1.2f TeV ",hd->GetMean());
+ text0.DrawLatex(xt,20.,ch);
+ char chd[60];
+ if(strncmp(system,"ZNA",3) == 0) sprintf(chd,"<E_{ZNA}> Glauber = %1.2f TeV",hg->GetMean());
+ else sprintf(chd,"<E_{ZPA}> 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);
+
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff