Updating
authorcoppedis <Chiara.Oppedisano@to.infn.it>
Thu, 17 Jul 2014 10:00:39 +0000 (12:00 +0200)
committercoppedis <Chiara.Oppedisano@to.infn.it>
Thu, 17 Jul 2014 10:00:39 +0000 (12:00 +0200)
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h
PWGPP/EVCHAR/GlauberSNM/makeCentralityFitnozeri.C

index 3e003a4..c6306e8 100755 (executable)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*   Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch                   */
-
-///////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-//  class to determine centrality percentiles from 1D distributions          // 
-//                                                                           //
-///////////////////////////////////////////////////////////////////////////////
-
 #include <TNamed.h>
 #include <TH1D.h>
 #include <TString.h>
 #include <TFile.h>
 #include <TMath.h>
+#include <TRandom1.h>
 #include <TROOT.h>
 #include <TH2F.h>
 #include <TF1.h>
@@ -325,8 +318,8 @@ void AliCentralityGlauberFit::MakeFits()
 
 
 //--------------------------------------------------------------------------------------------------
-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) 
+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");
@@ -335,10 +328,12 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
   
   TFile *outFile = NULL;
   TNtuple *ntuple = NULL;
+  
+  Double_t ntot=0, nblack=0, ngray=0, Etot=0, lcp=0, nslowp=0, nslown=0;
+  
   if(save){
     outFile = new TFile(fOutntuplename,"RECREATE");
-    ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot:nblack:ngray:Etot:lcp:nslowp:nslown");
+    ntuple = new TNtuple("gnt", "Glauber ntuple", "fNpart:fNcoll:fB:fTaa:ntot:nblack:ngray:Etot:lcp:nslowp:nslown");
   } 
 
   Int_t nents=0, evtot=0, evn=0, evzeron=0, evzerop=0, evenzero=0;
@@ -351,24 +346,30 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
       fNpart = 2;
       fNcoll = 1;
     }
-    //printf(" Getting entry %d from ntuple\n",i%nents);
+    //printf(" Getting entry %d from ntuple  Ncoll = %f\n",i%nents, fNcoll);
 
     // Slow Nucleon Model from Chiara
-    Double_t ntot=0., n=0., lcp=0., np=0., nn=0.;
-    Double_t nbn=0., ngn=0., nbp=0., ngp=0.;
+    Double_t n=0., nbn=0., ngn=0., nbp=0., ngp=0.;
+    //
     //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+nbn); // ZNA
-    else      n = alpha*(ngp+nbp); // ZPA
-    
-    np = nbp+ngp;
-    nn = nbn+ngn;
+    if(fIsZN){
+      n = alpha*(ngn+nbn); 
+      nblack = nbn;
+      ngray = ngn;
+    }
+    else{
+      n = alpha*(ngp+nbp); // ZPA
+      nblack = nbp;
+      ngray = ngp;
+    }    
+    nslown = nbn+ngn;
+    nslowp = nbp+ngp;
 
     if((alpha*(ngn+nbn)<0.5)) evzeron++;
     if((alpha*(ngp+nbp)<0.5)) evzerop++;
@@ -383,22 +384,14 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
       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);
-    //----------------------------------------
-
+    //else n=0.;
+    
     // non-lineary effect -------------------------------
-    ntot = ntot + k*ntot*ntot;
+    //ntot = ntot + k*ntot*ntot;
 
     Double_t nFact = 4.*82./208.;
-
-    int ntote = (int) n;
-    Float_t Etot = nFact*ntote;
+    Etot = nFact*n;
         
-    //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.){
@@ -406,19 +399,14 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
       //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;
+    //Etot = Etot + k*Etot*Etot;
 
     if(Etot<0.5) evenzero++;
      
-    if(Etot>0) {
-      h1->Fill(Etot);     
-      if(save) ntuple->Fill(fNpart, fNcoll, fB, fTaa, (float) ntot, nbn, ngn, Etot, lcp, np, nn);
-    }
+    h1->Fill(Etot);    
+    if(save) ntuple->Fill(fNpart, fNcoll, fB, fTaa, ntot, nblack, ngray, Etot, lcp, nslowp, nslown);
   }
 
   if(save) {
@@ -557,35 +545,37 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
   Float_t nu = (Float_t) ncoll; 
   //
   nu = gRandom->Gaus(nu,0.5); 
-  //if(nu<1.) nu = 1.;
+  if(nu<1.) nu = 1.;
 
+  // PROTONS ----------------------------------------------------------------
   //  gray protons
   Float_t  poverpd = 0.843; 
   Float_t  zAu2zPb = 82./79.;
   Float_t  nGrayp = (-0.27 + 0.63 * nu - 0.0008 *nu *nu)*poverpd*zAu2zPb;
+  if(nGrayp<0.) nGrayp=0;
 
   Double_t p;
   p =  nGrayp/fP;
-  ngp = gRandom->Binomial((Int_t) fP, p);
-  //if(nGrayp<0.) ngp=0;
+  ngp = (double) gRandom->Binomial(fP, p);
 
   //  black protons
   //Float_t blackovergray = 3./7.;// 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 = gRandom->Binomial((Int_t) fP, p);
-  //
-  //if(nBlackp<0.) nbp=0;
+  p =  nBlackp/(fP-ngp);
+  nbp = (double) gRandom->Binomial((Int_t) (fP-ngp), p);
+  // SATURATION in no. of black protons
+  //if(ngp>=9.) nbp = 15;
   
  // NEUTRONS ----------------------------------------------------------------
 
-  if(nu<3.){
+  /*if(nu<3.){
     nGrayp = -0.836 + 0.9112 *nu - 0.05381 *nu *nu;
     nBlackp  = blackovergray*nGrayp; 
-  }  
+  } */ 
 
   //  gray neutrons
   Float_t nGrayNeutrons = 0.;
@@ -594,20 +584,24 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
 
   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;
+//changed
+//    float xconj = fParamnSlow[1]/fParamnSlow[0]-fParamnSlow[2];
+    float xconj = 1.;
+    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;  // now integrated in fParamnSlow[0], fParamnSlow[1]
     //
-    nBlackNeutrons = nSlow * 0.9;;
+    nBlackNeutrons = nSlow * 0.9;
+    //nBlackNeutrons = nSlow * 0.5;
     nGrayNeutrons  = nSlow - nBlackNeutrons;
   }
   else{
     // slightly adjusted Sikler corrected for neutron yield...
-    nBlackNeutrons = 126./208. * 3.6 * nu; 
+    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>)
@@ -617,18 +611,26 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
   if(nBlackNeutrons<0.) nBlackNeutrons=0.;
 
   p =  nGrayNeutrons/fN;
-  ngn = gRandom->Binomial((Int_t) fN, p);
+  ngn = gRandom->Binomial(fN, 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)));
+  p =  nBlackNeutrons/(fN-ngn);
+  nbn = gRandom->Binomial((Int_t) (fN-ngn), 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.;
   
 }
 
 //--------------------------------------------------------------------------------------------------
- void AliCentralityGlauberFit::MakeSlowNucleons2s(Int_t ncoll, Double_t bog, Double_t gamma, Double_t sigmap, 
+ void AliCentralityGlauberFit::MakeSlowNucleons2s(Float_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)
@@ -636,10 +638,11 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
 // Return the number of black and gray nucleons
 //
 
-  Int_t fP = 82;
-  Int_t fN = 126;
+  float fP = 82;
+  float fN = 126;
   
-  Float_t nu = (Float_t) ncoll; 
+  Float_t nu = ncoll;
+  //Float_t nu = ncoll*(gRandom->Rndm()+0.5);
 
  // PROTONS ----------------------------------------------------------------
   //  gray protons
@@ -648,6 +651,8 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
   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;
   //
+//changed
+  //Float_t  nGrayp = grayp;
   Float_t  nGrayp = gRandom->Gaus(grayp, sigmap);
 //changed
 //  Float_t  nGrayp = gRandom->Gaus(grayp, (0.8*sigmap-sigmap)*grayp/90+sigmap);
@@ -655,7 +660,7 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
 
   Double_t p=0.;
   p =  nGrayp/fP;
-  ngp = (double) gRandom->Binomial(fP, p);
+  ngp = (double) gRandom->Binomial((int) fP, p);
 
   //  black protons
   //Float_t blackovergray = 3./7.;// =0.43 from spallation
@@ -663,18 +668,21 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
   Float_t blackovergray = bog;
   //
   Float_t nBlackp = blackovergray*nGrayp;
-  if(nBlackp<0.) nBlackp=0;
+  //if(nBlackp<0.) nBlackp=0;
 
-  p =  nBlackp/fP;
-  nbp = (double) gRandom->Binomial(fP, p);
+  p =  nBlackp/(fP-ngp);
+  nbp = (double) gRandom->Binomial((int) (fP-ngp), p);
+  // SATURATION in no. of black protons
+  //if(ngp>=9.) nbp = 15;
 
  // NEUTRONS ----------------------------------------------------------------
   // Relevant ONLY for n
-  // NB DOESN'T WORK FOR SMEARING IN Nslow!!!!!!!!!!
+  // 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);
+    //nGrayp = gRandom->Gaus(grayp, sigmap);
+    nGrayp = grayp;
     nBlackp  = blackovergray*nGrayp; 
   }*/
     
@@ -689,14 +697,15 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
     //
 //changed
 //    float xconj = fParamnSlow[1]/fParamnSlow[0]-fParamnSlow[2];
-    float xconj = 2.;
+    /*float xconj = 1.;
     float yconj = fParamnSlow[0]-fParamnSlow[1]/(fParamnSlow[2]+xconj)+fParamnSlow[3]*xconj;
-    if(lcp<xconj) nSlow = 0.+(yconj-0.)*(lcp-0.)/(xconj-0.);
+    if(lcp<xconj) nSlow = 0.+(yconj-0.)*(lcp-0.)/(xconj-0.);*/
     //
     // Factor to scale COSY to obtain ALICE n mltiplicities
-    nSlow = 1.68*nSlow;
+    //    nSlow = 1.68*nSlow;  // now integrated in fParamnSlow[0], fParamnSlow[1]
     //
-    nBlackNeutrons = nSlow * 0.9;;
+    nBlackNeutrons = nSlow * 0.9;
+    //nBlackNeutrons = nSlow * 0.5;
     nGrayNeutrons  = nSlow - nBlackNeutrons;
   }
   else{
@@ -721,126 +730,28 @@ void AliCentralityGlauberFit::MakeSlowNucleons2(Int_t ncoll, Double_t bog, Doubl
   if(nGrayNeutrons<0.) nGrayNeutrons=0.;
   if(nBlackNeutrons<0.) nBlackNeutrons=0.;
 
-  p =  nGrayNeutrons/fN;
-  //ngn = gRandom->Binomial((Int_t) fN, 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->Binomial((Int_t) fN, p);
+  TRandom1 r;
+  //nbn = (double) r.Binomial((int) fN, 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(nBlackNeutrons>=10) nbn = r.Gaus(nBlackNeutrons+0.5, TMath::Sqrt(fN*p*(1-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::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)
-//
-// Return the number of black and gray nucleons
-//
-// Number of collisions
-
-   // based on E910 model ================================================================
-
-  Int_t fP = 82;
-  Int_t fN = 126;
-  
-  Float_t nu = (Float_t) ncoll; 
-  Float_t sigmap = 0.25;
-
-  //  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;
-  Float_t  nGrayp = gRandom->Gaus(grayp, sigmap);
-  //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)));
-
-  //  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);
-  Float_t nBlackp = blackovergray*nGrayp;
-  //if(nBlackp<0.) nBlackp=0.;
-
-  p =  nBlackp/fP;
-  nbp = gRandom->Binomial((Int_t) fP, p);
-  
   //  gray neutrons
-  Float_t nGrayNeutrons = 0.;
-  Float_t nBlackNeutrons = 0.;
-  Float_t cp = (nGrayp+nBlackp)/gamma;
+  p =  nGrayNeutrons/(fN-nbn);
+  ngn = gRandom->Binomial((Int_t) (fN-nbn), p);
+  //ngn = (double) (r.Binomial((Int_t) (fN-nbn), p));
+//changed
+  /*if(nGrayNeutrons>=10) ngn = gRandom->Gaus(nGrayNeutrons+0.5, TMath::Sqrt(fN*p*(1-p)));
+  else ngn = gRandom->Binomial((Int_t) fN, p);*/
+  if(ngn<0.) ngn=0.;
 
-  if(cp>0.){
-    //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;
-    nBlackNeutrons = nSlow - nGrayNeutrons;
-  }
-  else{
-    // Sikler "pasturato"
-    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);
-  }
-  //
-  if(nGrayNeutrons<0.) nGrayNeutrons=0.;
-  if(nBlackNeutrons<0.) nBlackNeutrons=0.;
+  //printf(" Ncoll %1.2f   Ngp %1.2f Nbp %1.2f   Nbn %1.2f Ngn %1.2f\n",nu, ngp,nbp,nbn,ngn);
 
-  p =  nGrayNeutrons/fN;
-  //ngn = gRandom->Binomial((Int_t) fN, p);
-  ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN*p*(1-p)));
+}
 
-  //  black neutrons
-  p =  nBlackNeutrons/fN;
-  //nbn = gRandom->Binomial((Int_t) fN, p);
-  nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p)));
-  
-   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)
 { 
index dcda9f9..9b0e21b 100755 (executable)
@@ -49,42 +49,42 @@ class AliCentralityGlauberFit : public TObject {
   void     SetChi2Min(double chi2min) {fChi2Min=chi2min;}
   
  private:
-  Int_t   fNk;            // k points
+  Int_t    fNk;           // k points
   Double_t fKlow;         // k low
   Double_t fKhigh;        // k high
-  Int_t   fNalpha;        // alpha points
+  Int_t    fNalpha;       // alpha points
   Double_t fAlphalow;     // alpha low
   Double_t fAlphahigh;    // alpha high
-  Int_t   fNsigma;        // sigma points
+  Int_t    fNsigma;       // sigma points
   Double_t fSigmalow;     // sigma low
   Double_t fSigmahigh;    // sigma high
-  Int_t   fNbog;          // bog points
+  Int_t    fNbog;         // bog points
   Double_t fBoglow;       // bog low
   Double_t fBoghigh;      // bog high
-  Int_t   fNgamma;        // gamma points
+  Int_t    fNgamma;       // gamma points
   Double_t fgammalow;     // gamma low
   Double_t fgammahigh;    // gamma high
-  Int_t   fNsigmap;       // sigmap points
+  Int_t    fNsigmap;      // sigmap points
   Double_t fsigmaplow;    // sigmap low
   Double_t fsigmaphigh;   // sigmap high
-  Int_t   fRebinFactor;   // rebin factor
+  Int_t    fRebinFactor;  // rebin factor
   Double_t fScalemin;     // mult min where to scale
   Double_t fMultmin;      // mult min
   Double_t fMultmax;      // mult max
   TNtuple *fGlauntuple;   // glauber ntuple
-  Float_t fNpart;         // number participants
-  Float_t fNcoll;         // number collisions
-  Float_t fB;             // impact parameter
-  Float_t fTaa;           // taa
-  TH1F  *fTempHist;       // Temporary pointer to data histo, to be used in minimization 
-  TH1F  *fGlauberHist;    // Glauber histogram
-  Bool_t fUseChi2;        // If true, use chi2
-  Int_t fNevents;         // Number of events to use in the glauber ntuple
-  Bool_t fApplySaturation;// 
-  Int_t fnGraySaturation; //
-  Int_t fnBlackSaturation;//
-  Bool_t  fIsZN;
-  Float_t fSDprobability[40];
+  Float_t  fNpart;        // number participants
+  Float_t  fNcoll;        // number collisions
+  Float_t  fB;            // impact parameter
+  Float_t  fTaa;          // taa
+  TH1F    *fTempHist;     // Temporary pointer to data histo, to be used in minimization 
+  TH1F    *fGlauberHist;     // Glauber histogram
+  Bool_t   fUseChi2;         // If true, use chi2
+  Int_t    fNevents;         // Number of events to use in the glauber ntuple
+  Bool_t   fApplySaturation; 
+  Int_t    fnGraySaturation; 
+  Int_t    fnBlackSaturation;
+  Bool_t   fIsZN;
+  Float_t  fSDprobability[40];
   Double_t fParamnSlow[4];
   Int_t    fNTrials;
   Double_t fChi2Min;
@@ -103,10 +103,8 @@ class AliCentralityGlauberFit : public TObject {
   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,
+  void     MakeSlowNucleons2s(Float_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);
 
index b298082..99624d5 100644 (file)
@@ -27,7 +27,7 @@
 #endif
 
 void makeCentralityFitnozeri(const int nRun=195483, const char *system = "ZNA",
-       int ntrials=1, int Rebin=1, int Nevt=1.e7)
+       int ntrials=1, int Rebin=1, int Nevt=1.e6)
 {
  //load libraries
   gSystem->SetBuildDir(".");
@@ -80,9 +80,8 @@ void makeCentralityFitnozeri(const int nRun=195483, const char *system = "ZNA",
   // 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);
+  mPM->SetNParam(50., 230., 4.4, 0.48);
+  //mPM->SetNParam(40., 230., 4.5, 0.45);
   //
   // ALICE ex-novo
   //mPM->SetNParam(78., 700., 7.);
@@ -91,18 +90,22 @@ void makeCentralityFitnozeri(const int nRun=195483, const char *system = "ZNA",
   if(strncmp(system,"ZNA",3) == 0) {
     printf(" Setting parameters for ZNA Glauber fit\n\n");
     mPM->SetIsZN();
-    mPM->SetRangeToFit(0.5, 90.5);   
+    mPM->SetRangeToFit(0.8, 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);
+    // original
+    if(nRun==195483) mPM->SetGlauberParam(1,0.0,0.0, 1,0.957,1., 1,0.26,0.30, 1,0.65,0.65, 1,0.585,0.585, 1, 0.25,0.3);
+    //to study systematics
+    //if(nRun==195483) mPM->SetGlauberParam(1,0.0,0.0, 1,0.956,1., 1,0.29,0.30, 1,1.80,0.65, 1,0.585,0.585, 1, 0.25,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.585,0.585, 1, 0.25,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);
+    // original
+    //if(nRun==195483) mPM->SetGlauberParam(1,0.0,0.0, 1,0.60,0.8, 1,0.6,0.65, 1,0.65,0.65, 1,0.585,0.585, 1, 0.25,0.4);
+    if(nRun==195483) mPM->SetGlauberParam(1,0.0,0.0, 1,0.40,0.8, 1,0.4,0.65, 1,1.80,0.65, 1,0.585,0.585, 1, 0.25,0.4);
+    else             mPM->SetGlauberParam(1,0.0,0.0, 1,0.60,0.8, 1,0.5,0.65, 1,0.65,0.65, 1,0.585,0.585, 1, 0.25,0.4);
   }
   mPM->MakeFits();  
 
@@ -187,6 +190,7 @@ void makeCentralityFitnozeri(const int nRun=195483, const char *system = "ZNA",
   if(strncmp(system,"ZNA",3) == 0) sprintf(psn,"ZNA_fit%d.gif",nRun);
   else  sprintf(psn,"ZPA_fit%d.gif",nRun);
   g->Print(psn);
-
+       
+       printf("\n Everything is OK: ntuple and fit results saved!!! \n\n");
 }