Updated SNM Glauber fit
authorcoppedis <chiara.oppedisano@cern.ch>
Thu, 3 Jul 2014 10:59:14 +0000 (12:59 +0200)
committercoppedis <chiara.oppedisano@cern.ch>
Thu, 3 Jul 2014 11:06:37 +0000 (13:06 +0200)
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h
PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C [new file with mode: 0644]

index a95f235..3e003a4 100755 (executable)
@@ -34,7 +34,7 @@
 #include <TGraphErrors.h>
 #include <vector>
 #include <TMinuit.h>
-#include <TStopwatch.h>
+#include <TCanvas.h>
 #include <TRandom.h>
 #include <TDatabasePDG.h>
 #include <TPDGCode.h>
@@ -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<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
@@ -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;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;
 }
 
@@ -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 = 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;
@@ -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<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)
 //
@@ -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)));
     
index 31fafe6..dcda9f9 100755 (executable)
@@ -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<TString> 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 (file)
index 0000000..b298082
--- /dev/null
@@ -0,0 +1,192 @@
+#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);
+
+}
+