removed smear in ntot

diff --git a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
index 7e674a3fe0cdd87dde7d1a7b2db5afad06ec7701..a95f2359560e2a136f6a07cb57c41a110408192f 100755 (executable)
--- a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
+++ b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
@@ -78,7 +78,8 @@ AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) :
   fOutrootfilename(0),
   fOutntuplename(0),
   fAncfilename("ancestor_hists.root"),
-  fHistnames()
+  fHistnames(),
+  fIsZN(kTRUE)
 {
   // Standard constructor.
   TFile *f = 0;
@@ -167,12 +168,13 @@ 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;
     if (!hDATA) {
       TList *list  = (TList*) (inrootfile->Get("CentralityStat")); 
       //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat")); 
       hDATA  = (TH1F*) (list->FindObject(*hni));
     } 
-    hDATA->Rebin(fRebinFactor);
+    //hDATA->Rebin(fRebinFactor);
     //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
 
     Double_t chi2min = 9999999.0;
@@ -214,6 +216,7 @@ void AliCentralityGlauberFit::MakeFits()
        }
       }
     }
+    thistGlau->Reset();
     thistGlau = GlauberHisto(k_min,alpha_min,sigma_min,bog_min,CP_min, hDATA,kTRUE);
 
     TH1F * hGLAU = 0x0;
@@ -222,13 +225,27 @@ void AliCentralityGlauberFit::MakeFits()
     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f_%.3f_%.3f",
                                                              k_min,alpha_min,sigma_min,bog_min,CP_min)));
 
-    Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
-    Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
-    std::cout << hGLAU->FindBin(fScalemin) << " " << hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX()) << std::endl;
-    std::cout << hDATA->FindBin(fScalemin) << " " << hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX()) << std::endl;
+    Int_t lastBin=0;
+    if(fIsZN) lastBin = 100;//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;
+    //
+    //printf(" binwidth: hGLAU %f  hDATA %f\n", hGLAU->GetBinWidth(100), hDATA->GetBinWidth(100));
+    printf(" scale %f \n", scale);
+    
     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;
 
-    SaveHisto(hDATA,hGLAU,outrootfile);
+    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);
@@ -251,7 +268,7 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
                                             TH1F *hDATA, Bool_t save) 
 {
   // Get Glauber histogram.
-  static TH1F *h1 = (TH1F*)hDATA->Clone();
+  static TH1F *h1 = (TH1F*)hDATA->Clone("h1");
   h1->Reset();
   h1->SetName(Form("fit_%.3f_%.3f",k,alpha));
 
@@ -260,7 +277,7 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
  
   if (save) {
     outFile = new TFile(fOutntuplename,"RECREATE");
-    ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nbn:ngn");
+    ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nbn:ngn:Etot");
   } 
 
   Int_t nents = 0;
@@ -276,36 +293,15 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
     }
 
     // Slow Nucleon Model from Chiara
-    Double_t n=0.;
-    Double_t ntot=0., 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);
-
-
-    //conversion to energy *13.5; // conversion to energy------------------------
-    //Double_t  Engn, Enbn;
-    //Engn = ConvertToEnergy(0.05);
-    //Enbn = ConvertToEnergy(0.005);
-
-    //Engn = 1;//gRandom->Gaus(13.5, TMath::Sqrt(13.5));
-    //Enbn = 1;//gRandom->Gaus(13.5, TMath::Sqrt(13.5));
-
-    //Enbn =(int)nbn*Enbn;
-    //Engn =(int)ngn*Engn;
-
-    //Enbn=0;
-    //for (int i=0; i<(int)nbn; i++){
-    // Enbn += ConvertToEnergy(0.005);
-    //}
-    //for (int i=0; i<(int)ngn; i++){
-    //  Engn += ConvertToEnergy(0.05);
-    // }
-
-    //cout << "GRAY " << Engn << " BLACK " << Enbn << endl;
+    Double_t ntot=0., n=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
-    n=alpha*ngn+alpha*nbn; // ZNA
-    //n=alpha*ngp+alpha*nbp; // ZPA
+    //
+    if(fIsZN) n = alpha*ngn+alpha*nbn; // ZNA
+    else n =  alpha*ngp+alpha*nbp; // ZPA
     //----------------------------------------
 
     //------ experimental resolution -> Gaussian smearing ------------------------------------
@@ -313,22 +309,32 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
     //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_t)(gRandom->Gaus(n, resexp));
+        
+    //if(n>0) resexp = sigma*TMath::Sqrt(n);    
+    //else resexp=0;
+    ntot = (int) (gRandom->Gaus(n, resexp));
     //----------------------------------------
 
     // non-lineary effect -------------------------------
     //ntot = k*ntot;
-    ntot = ntot + k*ntot*ntot;
+    //ntot = ntot + k*ntot*ntot;
 
     //cout << ngn << " " << nbn << " "  << ntot << endl;
 
-    if (ntot>0) {
-      h1->Fill(ntot);     
+    Double_t nFact = 1.577;
+    Double_t Etot = nFact*ntot;
+    //
+
+    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) {
+      h1->Fill(Etot);     
       if (save) 
-       ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot,nbn,ngn);
+       ntuple->Fill(fNpart, fNcoll, fB, fTaa, ntot, nbn, ngn, Etot);
     }
   }
 
@@ -348,8 +354,8 @@ Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau)
   Int_t lowchibin =   hDATA->FindBin(fMultmin);
   Int_t highchibin =  hDATA->FindBin(fMultmax);
 
-  Double_t mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
-  Double_t scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral);
+  Double_t mcintegral = thistGlau->Integral(1, thistGlau->GetNbinsX());
+  Double_t scale = (hDATA->Integral(1, hDATA->GetNbinsX())/mcintegral);
   thistGlau->Scale(scale);
 
   // calculate the chi2 between MC and real data over some range ????
@@ -464,6 +470,8 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   nu = gRandom->Gaus(nu,0.5); 
   if(nu<1.) nu = 1.;
   //
+  //float sdp = gRandom->Rndm();
+  //if(nu==1. && sdp<=0.2) return;
 
   //  gray protons
   Float_t  poverpd = 0.843; 

diff --git a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h
index e10e2b42301b6c55971dbb6fc10a3b72127c7402..31fafe6f8fe55a7c69daa55aecad07cbe32f5709 100755 (executable)
--- a/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h
+++ b/PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h
@@ -40,7 +40,9 @@ class AliCentralityGlauberFit : public TObject {
   void     UseChi2(Bool_t f)  { fUseChi2=f; }
   void     SetSaturation(Bool_t saturation) {fApplySaturation = saturation;}
   void     SetSaturationParams(Int_t ngray=15, Int_t nblack=28) {fnGraySaturation=ngray; fnBlackSaturation=nblack;}
-
+  void     SetIsZN() {fIsZN=kTRUE;}
+  void     SetIsZP() {fIsZN=kFALSE;}
+  
  private:
   Int_t   fNk;            // k points
   Double_t fKlow;         // k low
@@ -73,6 +75,7 @@ class AliCentralityGlauberFit : public TObject {
   Bool_t fApplySaturation;// 
   Int_t fnGraySaturation; //
   Int_t fnBlackSaturation;//
+  Bool_t fIsZN;
 
   TString fInrootfilename;          // input root file
   TString fInntuplename;            // input Gauber ntuple
@@ -94,6 +97,6 @@ class AliCentralityGlauberFit : public TObject {
   AliCentralityGlauberFit(const AliCentralityGlauberFit&);
   AliCentralityGlauberFit &operator=(const AliCentralityGlauberFit&);
 
-  ClassDef(AliCentralityGlauberFit, 1)  
+  ClassDef(AliCentralityGlauberFit, 2)  
 };
 #endif

diff --git a/PWGPP/EVCHAR/GlauberSNM/makeCentralityFit.C b/PWGPP/EVCHAR/GlauberSNM/makeCentralityFit.C
index 62a0b7f5e96518e57de59fa3770cc98232503ac1..1a8c18990e9335bd595f31a0446e9ed8a9b8da90 100755 (executable)
--- a/PWGPP/EVCHAR/GlauberSNM/makeCentralityFit.C
+++ b/PWGPP/EVCHAR/GlauberSNM/makeCentralityFit.C
@@ -1,4 +1,32 @@
-void makeCentralityFit(const char * run="188359",const char * system="ZNA", int Rebin=1,int Nevt=1e6)
+#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 makeCentralityFit(const char * run="188359", const char *system = "ZNA", int Rebin=1, int Nevt=1e6)
 {
  //load libraries
   gSystem->SetBuildDir("/tmp/");
@@ -12,16 +40,26 @@ void makeCentralityFit(const char * run="188359",const char * system="ZNA", int
   gROOT->LoadMacro("AliCentralityGlauberFit.cxx+");
 
   const char *finnameGlau ="GlauberMC_pPb_ntuple_sigma70_mind4_r662_a546_Rpro6.root";
-
-  // data
+  char histname[8];
+  
+  if((strncmp(system,"ZNA",3))==0){
+    printf("\n  Glauber fit on ZNA spectrum\n\n");
+    sprintf(histname,"hZNA");
+  }
+  else if((strncmp (system,"ZNs",3)) == 0){
+    printf("\n  Glauber fit on ZNA spectrum subtracting ZNC contribution (~SD)\n\n");
+    sprintf(histname,"hZPA");
+  }
+  else if((strncmp (system,"ZPA",3)) == 0){
+    printf("\n  Glauber fit on ZPA spectrum\n\n");
+    sprintf(histname,"hZPA");
+  }
+  
   TString finname = Form("Histos%s.root",run);
   TString foutname = Form("%s_fit_%s.root",system,run);
   TString foutnameGlau = Form("%s_ntuple_%s.root",system,run);
 
-  const char *histname=("hnSpecn");
-  //const char *histname=("hnSpecp");
   
-
   AliCentralityGlauberFit *mPM = new AliCentralityGlauberFit(finnameGlau);
   mPM->SetInputFile(finname);        
   mPM->SetInputNtuple(finnameGlau);     
@@ -33,26 +71,64 @@ void makeCentralityFit(const char * run="188359",const char * system="ZNA", int
   mPM->SetNevents(Nevt);
   mPM->UseChi2(kTRUE);     // If TRUE minimize Chi2
 
-  if (strncmp (system,"ZNA",3) == 0) {
-    mPM->SetRangeToFit(1., 70.);   
-    mPM->SetRangeToScale(1.);  
-    mPM->SetGlauberParam(1,0.00,1., 1,0.97,1., 1,0.95,0.3, 1,0.65,0.8, 1,0.585,0.6);
+  if (strncmp(system,"ZNA",3) == 0) {
+    mPM->SetIsZN();
+    mPM->SetRangeToFit(1., 300.);   
+    mPM->SetRangeToScale(0);  
+    mPM->SetGlauberParam(1,0.,1., 1,0.96,1., 2,0.22,0.25, 1,0.65,0.7, 1,0.56,0.585);
   }
-  else if (strncmp (system,"ZPA",3) == 0) {
-    mPM->SetRangeToFit(1., 25.);  2
+  else if (strncmp(system,"ZPA",3) == 0) {
+    mPM->SetIsZP();
+    mPM->SetRangeToFit(1., 25.);  
     mPM->SetRangeToScale(1.);  
-    mPM->SetGlauberParam(1,0.0,1., 1,0.6,1, 1,0.25,0.3, 1,0.65,0.8, 1,0.585,0.6);
+    mPM->SetGlauberParam(1,0.0,1., 1,0.6,1, 1,0.25,0.3, 1,0.65,0.8, 1,0.58,0.59);
   }
   mPM->MakeFits();  
 
-  TFile * f = new TFile (foutname);
-  TH1 * hd = (TH1*) gDirectory->Get(("hnSpecn"));
-  TH1 * hg = (TH1*) gDirectory->Get(("hnSpecn_GLAU"));
-  hg->SetLineColor(kRed);
-  hd->Draw("hist");
-  hg->Draw("histsame");
-
+  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.);
+  hg->Draw("E");
+  hd->Draw("PEsame");
+  hd->SetXTitle("E_{ZNA} (TeV)");
+  hg->Draw("Esame");
+  gPad->SetLogy(1);
+  
+  TLatex text0;
+  text0.SetTextSize(0.04);
+  text0.SetTextColor(kBlue+3);
+  char ch[60];
+  sprintf(ch,"<E_{ZNA}> DATA = %1.2f TeV ",hd->GetMean());
+  text0.DrawLatex(xt,10.,ch);
+  char chd[60];
+  sprintf(chd,"<E_{ZNA}> Glauber = %1.2f TeV",hg->GetMean());
+  text0.SetTextColor(kPink-2);
+  text0.DrawLatex(xt,5.,chd);
 
+  char ct[60];
+  sprintf(ct,"RUN %s",run);
+  text0.SetTextColor(kAzure);
+  text0.DrawLatex(xt,1.,ct);
 
 }