smearing gray and black p,not ncoll
authoratoia <atoia@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 24 Sep 2013 21:55:57 +0000 (21:55 +0000)
committeratoia <atoia@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 24 Sep 2013 21:55:57 +0000 (21:55 +0000)
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.cxx
PWGPP/EVCHAR/GlauberSNM/AliCentralityGlauberFit.h

index 943a33f..7e674a3 100755 (executable)
@@ -231,6 +231,7 @@ void AliCentralityGlauberFit::MakeFits()
     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;
     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
@@ -275,9 +276,10 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
     }
 
     // Slow Nucleon Model from Chiara
-    Double_t n;
-    Double_t ntot, nbn, ngn, nbp, ngp;
-    MakeSlowNucleons2(fNcoll,alpha,k,bog,CP, nbn,ngn,nbp,ngp);
+    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------------------------
@@ -302,19 +304,19 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t k, Double_t alpha, Double_t
     //cout << "GRAY " << Engn << " BLACK " << Enbn << endl;
 
     // acceptance
-    n=alpha*ngn+alpha*nbn;
-    // n=alpha*ngp+alpha*nbp;
+    n=alpha*ngn+alpha*nbn; // ZNA
+    //n=alpha*ngp+alpha*nbp; // ZPA
     //----------------------------------------
 
     //------ experimental resolution -> Gaussian smearing ------------------------------------
     Double_t resexp=0;
-    //if (n>0) resexp = 1./TMath::Sqrt(n)*sigma*gRandom->Gaus(0.0,1.0);
+    //if (n>0) resexp = sigma*gRandom->Gaus(0.0,1.0)/TMath::Sqrt(n);
     //else resexp=0;
-    //ntot=n*(1+resexp);
+    //ntot = n*(1+resexp);
 
-    if (n>0)resexp = sigma*TMath::Sqrt(n)/2;    
+    if(n>0) resexp = sigma*TMath::Sqrt(n);    
     else resexp=0;
-    ntot = (Int_t)(gRandom->Gaus(n,resexp));
+    ntot = (Int_t)(gRandom->Gaus(n, resexp));
     //----------------------------------------
 
     // non-lineary effect -------------------------------
@@ -443,10 +445,8 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
 }
 
 
-
 //--------------------------------------------------------------------------------------------------
  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 alpha, Double_t k, Double_t bog, Double_t CP, Double_t &nbn, Double_t &ngn)
 {
 // from AliGenSlowNucleonModelExp (Chiara Oppedisano)
 //
@@ -461,7 +461,8 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   
   Float_t nu = (Float_t) ncoll; 
   //
-  nu = nu+1.*gRandom->Rndm();
+  nu = gRandom->Gaus(nu,0.5); 
+  if(nu<1.) nu = 1.;
   //
 
   //  gray protons
@@ -479,6 +480,7 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   //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;
@@ -495,7 +497,6 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   //  gray neutrons
   Float_t nGrayNeutrons = 0.;
   Float_t nBlackNeutrons = 0.;
-  //Float_t cp = (nGrayp+nBlackp)/0.24;
   Float_t cp = (nGrayp+nBlackp)/CP;
 
   // if(cp>0.){
@@ -508,7 +509,7 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   //   nBlackNeutrons = nSlow - nGrayNeutrons;
   // }
   if(cp>0.){
-    Float_t paramnSlow[3] = {60., 469.2, 8.762};
+    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.);
@@ -518,9 +519,11 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   }
   else{
     // Sikler "pasturato"
-    nGrayNeutrons = 0.47 * 2.3 *  nu; // fAlphaGray=2.3
-    nBlackNeutrons = 0.88 * 3.6 * nu; // fAlphaBlack=3.6      
-    printf("nslowp=0 -> ncoll = %1.0f -> ngrayn = %1.0f  nblackn = %1.0f \n", nu, nGrayNeutrons, nBlackNeutrons);
+    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);
   }
 
   p =  nGrayNeutrons/fN;
@@ -533,6 +536,85 @@ void AliCentralityGlauberFit::MakeSlowNucleons(Int_t ncoll, Double_t alpha, Doub
   nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN*p*(1-p)));
 }
 
+//--------------------------------------------------------------------------------------------------
+ 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)
+{
+// 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)));
+  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.;
+
+  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;
+
+  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.);
+    
+    nGrayNeutrons = nSlow * 0.1;
+    nBlackNeutrons = nSlow - nGrayNeutrons;
+  }
+  else{
+    // Sikler "pasturato"
+    nGrayNeutrons = 0.47 * 2.2 *  nu; // fAlphaGray=2.3
+    nBlackNeutrons = 0.88 * 3.6 * nu; // fAlphaBlack=3.6  
+    //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.;
+
+  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;
+}
+
 Double_t AliCentralityGlauberFit::ConvertToEnergy(Double_t T)
 { 
    TDatabasePDG * pdg = TDatabasePDG::Instance();
index b987040..e10e2b4 100755 (executable)
@@ -86,6 +86,7 @@ class AliCentralityGlauberFit : public TObject {
   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 ConvertToEnergy(Double_t T);
   Double_t Maxwell(Double_t m, Double_t p, Double_t T);