update for pA

diff --git a/PWGPP/EVCHAR/AliCentralityBy1D.cxx b/PWGPP/EVCHAR/AliCentralityBy1D.cxx
index e655a1477b145860b129f2590b573020cda0bb4e..ec6b618ab7ae0c15d4da2d8491b2a216fae98f91 100644 (file)
--- a/PWGPP/EVCHAR/AliCentralityBy1D.cxx
+++ b/PWGPP/EVCHAR/AliCentralityBy1D.cxx
@@ -80,6 +80,7 @@ void AliCentralityBy1D::MakePercentiles() {
      hpercent->Reset();
 
      int start_bin=htemp->FindBin(fMultLowBound);
+     //cout << "START: " << start_bin << " " << fMultLowBound << " " << fPercentXsec << endl;
      for (int ibin=1; ibin<=htemp->GetNbinsX(); ibin++) {
       
        if (ibin>=start_bin)
@@ -88,7 +89,9 @@ void AliCentralityBy1D::MakePercentiles() {
                                htemp->Integral(start_bin,htemp->GetNbinsX()));
        else
        hpercent->SetBinContent(ibin, 100);
-     }    
+ 
+       //if (ibin<700) cout << ibin << " " << hpercent->GetBinContent(ibin) << endl;
+    }    
 
      SaveHisto(htemp,hpercent,outrootfile);
 

diff --git a/PWGPP/EVCHAR/AliCentralityGlauberFit.cxx b/PWGPP/EVCHAR/AliCentralityGlauberFit.cxx
index c839df97607bcba69e226d66fff27316f868c2b0..771abf742ebbd104877b493ecf9978a86f4c9494 100644 (file)
--- a/PWGPP/EVCHAR/AliCentralityGlauberFit.cxx
+++ b/PWGPP/EVCHAR/AliCentralityGlauberFit.cxx
@@ -159,6 +159,7 @@ void AliCentralityGlauberFit::MakeFits()
     hDATA  = (TH1F*) (inrootfile->Get(*hni)); 
     if (!hDATA) {
       TList *list  = (TList*) (inrootfile->Get("CentralityStat")); 
+      //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat")); 
       hDATA  = (TH1F*) (list->FindObject(*hni));
     } 
     hDATA->Rebin(fRebinFactor);
@@ -249,7 +250,13 @@ void AliCentralityGlauberFit::MakeFitsMinuitNBD(Double_t alpha, Double_t mu, Dou
   std::vector<TString>::const_iterator hni;
   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
     hDATA  = (TH1F*) (inrootfile->Get(*hni)); 
+    if (!hDATA) {
+      TList *list  = (TList*) (inrootfile->Get("CentralityStat")); 
+      //TList *list  = (TList*) (inrootfile->Get("VZEROEquaFactorStat")); 
+      hDATA  = (TH1F*) (list->FindObject(*hni));
+    } 
     hDATA->Rebin(fRebinFactor);
+
     fTempHist=hDATA;
     TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
     hGLAU->Sumw2();
@@ -344,11 +351,11 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t mu, Double_t k, Double_t al
       Double_t weights = fhAncestor->GetBinContent(np);
 
       if (weights <= 0) continue;
-      Int_t trials = (Int_t) (20 * nanc * (int) mu);
+      Int_t trials = (Int_t) (20 * nanc * (Int_t) mu);
       if (trials <=0) continue;
       for (Int_t j=0; j<trials; j++) {
-               double nbdvalue = NBD(j, mu * nanc, k * nanc);
-               h1->Fill((double) j, nbdvalue * weights);
+               Double_t nbdvalue = NBD(j, mu * nanc, k * nanc);
+               h1->Fill((Double_t) j, nbdvalue * weights);
       }
     }
     return h1;
@@ -376,8 +383,8 @@ TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t mu, Double_t k, Double_t al
       fNcoll = 1;
     }
     Int_t n=0;
-    //if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNpart,alpha));
-    if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNcoll,alpha));
+    if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNpart,alpha));
+    //if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNcoll,alpha));
     else if (fAncestor == 2) n = (Int_t)(alpha * fNpart + (1-alpha) * fNcoll);
     else if (fAncestor == 3) n = (Int_t)((1-alpha) * fNpart/2 + alpha * fNcoll);
 
@@ -492,8 +499,8 @@ void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TH1F *heffi, T
 Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) const
 {
   // Compute NBD.
-  double F;
-  double f;
+  Double_t F;
+  Double_t f;
 
   if (n+k > 100.0) {
     // log method for handling large numbers
@@ -536,11 +543,11 @@ void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t
   Double_t mu    = par[1];
   Double_t k     = par[2];
 
-  if (0) { //avoid warning
-    gin=gin;
-    npar=npar;
-    iflag=iflag;
-  }
+  // if (0) { //avoid warning
+  //   gin=gin;
+  //   npar=npar;
+  //   iflag=iflag;
+  // }
   AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
   TH1F * thistGlau = obj->GlauberHisto(mu,k,alpha,obj->GetTempHist(),kFALSE);
   f = obj->CalculateChi2(obj->GetTempHist(),thistGlau);

diff --git a/PWGPP/EVCHAR/makeCentralityBy1D.C b/PWGPP/EVCHAR/makeCentralityBy1D.C
index c1d66cf4b9d13ec3f3598196e3131c814fa37362..e8970e7af7c35c08538a862bc97ad0841ab43a4e 100644 (file)
--- a/PWGPP/EVCHAR/makeCentralityBy1D.C
+++ b/PWGPP/EVCHAR/makeCentralityBy1D.C
@@ -1,7 +1,13 @@
 void makeCentralityBy1D 
-//(int run =167693, const char *histname ="fHOutMultV0M", float percentXsec=90.05, float mult=77.4) //V0M
-//(int run =167693, const char *histname ="fHOutMultCL1", float percentXsec=90.14, float mult=20) // CL1
-(int run =167693, const char *histname ="fHOutMultTRK", float percentXsec=90.09, float mult=10) // TRK
+//(int run =167987, const char *system ="V0M", float percentXsec=90.06, float mult=77.4) //V0M
+//(int run =170162, const char *system ="CL1", float percentXsec=90.04, float mult=20) // CL1
+//(int run =167693, const char *system ="TRK", float percentXsec=90.20, float mult=11) // TRK
+
+// for pA pilot run
+// for MC use percentXsec=100 and mult=0
+// for MC available also NPA
+// for pA use V0M, V0A, V0C, CL1, TRK, CND, (tbd FMD: not available in 146079) 
+(int run =195344, const char *system ="FMD", float percentXsec=100.0, float mult=0.0) 
 {
  //load libraries
   // gSystem->SetBuildDir("/tmp/");
@@ -15,8 +21,9 @@ void makeCentralityBy1D
   gROOT->LoadMacro("AliCentralityBy1D.cxx+");
   AliCentralityBy1D *mPM = new AliCentralityBy1D();
 
-  TString finname = Form("/home/atoia/analysis/data2011/multRef/AnalysisResults_%i.root",run);
-  TString foutname = Form("/home/atoia/analysis/data2011/by1D/AliCentralityBy1D_%i_TRK.root",run);
+  TString finname = Form("/home/atoia/analysis/data2013/multRef/EventStat_temp_%i.root",run);
+  TString foutname = Form("/home/atoia/analysis/data2013/by1D/AliCentralityBy1D_%i_%s.root",run,system);
+  const char *histname=Form("fHOutMult%s",system);
 
   mPM->AddHisto(histname);  
   mPM->SetInputFile(finname);        

diff --git a/PWGPP/EVCHAR/makeCentralityFit.C b/PWGPP/EVCHAR/makeCentralityFit.C
index 8de50d70b4a4f16d3c4effe7acc726ecb5e34619..1bfe72ce97fcb72e28d9956c994cc8caa9b10a6b 100644 (file)
--- a/PWGPP/EVCHAR/makeCentralityFit.C
+++ b/PWGPP/EVCHAR/makeCentralityFit.C
@@ -1,4 +1,4 @@
-void makeCentralityFit()
+void makeCentralityFit(const char * run="167693",const char * system="V0M", int Rebin=100,int Nevt=1e5, int Ncell=0)
 {
  //load libraries
   gSystem->SetBuildDir("/tmp/");
@@ -13,17 +13,19 @@ void makeCentralityFit()
 
   const char *finnameGlau ="/home/atoia/GlauberNtuple/GlauberMC_PbPb_ntuple_sigma64_mind4_r662_a546.root";
 
-  const char * run="167693";
-  const char * system="TRK";
-  TString finname = Form("/home/atoia/analysis/data2011/multRef/AnalysisResults_%s.root",run);
-  TString foutname = Form("/home/atoia/analysis/data2011/fit/%s_fit_%s.root",system,run);
-  //TString foutnameGlau = Form("/home/atoia/analysis/data2011/fit/%s_ntuple_%s.root",system,run);
-  TString foutnameGlau = Form("test_%s.root",system,run);
-  //const char *histname="fHOutMultV0M";
-  //const char *histname="fHOutMultCL1";
-  const char *histname="fHOutMultTRK";
 
+  TString finname = Form("/home/atoia/analysis/data2011/multRef/EventStat_temp_%s.root",run);
+  TString foutname = Form("/home/atoia/analysis/data2011/fit/%s_fitTEST_%s.root",system,run);
+  TString foutnameGlau = Form("/home/atoia/analysis/data2011/fit/%s_ntupleTEST_%s.root",system,run);
+  const char *histname=Form("fHOutMult%s",system);
 
+
+  /*
+  TString finname = Form("/home/atoia/analysis/EPVzero/VZEROEquaFactorStat.root");
+  TString foutname = Form("/home/atoia/analysis/data2011/fit/%sCell%d_fit_%s.root",system,Ncell,run);
+  TString foutnameGlau = Form("/home/atoia/analysis/data2011/fit/%sCell%d_ntuple_%s.root",system,Ncell,run);
+  const char *histname=Form("fMultCell_%d",Ncell);
+  */
   AliCentralityGlauberFit *mPM = new AliCentralityGlauberFit(finnameGlau);
   mPM->SetInputFile(finname);        
   mPM->SetInputNtuple(finnameGlau);     
@@ -31,35 +33,42 @@ void makeCentralityFit()
   mPM->SetOutputNtuple(foutnameGlau);
   mPM->AddHisto(histname);
 
-  mPM->SetRebin(1);
-  mPM->SetAncestorMode(2);   // If 1 use Npart**alpha, if 2 use alpha*Npart + (1-alpha)*Ncoll
-  mPM->SetFastFitMode(0);    // If 1 or 2 use cruder approximation to compute curve faster 1:NBD, 2:Gauss
-  mPM->UseChi2(kTRUE);       // If TRUE minimize Chi2
-  mPM->UseAverage(kFALSE);    // If TRUE use Average
+  mPM->SetRebin(Rebin);
+  mPM->SetNevents(Nevt);
+  mPM->SetAncestorMode(2); // 1: Npart**alpha, 2: alpha*Npart + (1-alpha)*Ncoll
+  mPM->SetFastFitMode(0);  // 1:NBD, 2:Gauss
+  mPM->UseChi2(kTRUE);     // If TRUE minimize Chi2
+  mPM->UseAverage(kFALSE); // If TRUE use Average
   mPM->SetNtrials(1);
-  mPM->SetNevents(1e5);
 
   // ----------range to fit--------------
-  //mPM->SetRangeToFit(100., 20000.);   // V0M
-  //mPM->SetRangeToFit(40., 5200.);   // CL1
-  mPM->SetRangeToFit(10., 2400.);   // TRK
-  // ----------range to scale--------------
-  //mPM->SetRangeToScale(100.); // V0M  
-  //mPM->SetRangeToScale(40.); // CL1  
-  mPM->SetRangeToScale(10.);  //TRK 
-  // ----------initial parameters--------------
-  //mPM->SetGlauberParam(10,29,31, 10,0.7,1.5, 10,0.85,0.87);
-  //mPM->MakeFitsMinuitNBD(0.801,30.,1.2);          // initial parameters
-  //mPM->SetGlauberParam(1,28,28.5, 1,1.291,1.9, 1,0.801,0.82); // V0M
-  //mPM->SetGlauberParam(1,7.25,7.5, 1,1.291,1.9, 1,0.801,0.82); // Cl1
-  mPM->SetGlauberParam(1,3.5,3.7, 1,1.291,1.5, 1,0.801,0.81); // TRK
-  // ----------done--------------
+  if (strncmp (system,"V0M",1) == 0) {
+    mPM->SetRangeToFit(100., 22000.);   // range to fit
+    mPM->SetRangeToScale(100.); // range to scale
+    mPM->SetGlauberParam(1,28.7,29., 1,1.601,1.5, 1,0.80,0.805); // fit parameters
+  }
+  else if (strncmp (system,"CL1",1) == 0) {
+    mPM->SetRangeToFit(40., 5400.);   
+    mPM->SetRangeToScale(40.); 
+    mPM->SetGlauberParam(1,7.13,7.3, 1,1.217,1.9, 1,0.802,0.815); 
+  }
+  else if (strncmp (system,"TRK",1) == 0) {
+    mPM->SetRangeToFit(10., 2600.);   
+    mPM->SetRangeToScale(10.);  
+    mPM->SetGlauberParam(1,3.9,4.2, 1,1.3,2.5, 1,0.801,0.81);
+  }
 
   mPM->MakeFits();  
-  
+
+  // ----------for Minuit--------------
+  //mPM->MakeFitsMinuitNBD(0.8,28.,1.29);          // initial parameters
+
+
   TFile * f = new TFile (foutname);
-  TH1 * hd = (TH1*) gDirectory->Get("fHOutMultTRK");
-  TH1 * hg = (TH1*) gDirectory->Get("fHOutMultTRK_GLAU");
+  TH1 * hd = (TH1*) gDirectory->Get(Form("fHOutMult%s",system));
+  TH1 * hg = (TH1*) gDirectory->Get(Form("fHOutMult%s_GLAU",system));
+  //TH1 * hd = (TH1*) gDirectory->Get(Form("fMultCell_%d",Ncell));
+  //TH1 * hg = (TH1*) gDirectory->Get(Form("fMultCell_%d_GLAU",Ncell));
   hg->SetLineColor(kRed);
   hd->Draw("e");
   hg->Draw("same");