Add code for glauber fits

Author: ddobrigk

Common/Tools/Multiplicity/README.md

@@ -0,0 +1,17 @@
+# Multiplicity/centrality tools in O2/O2Physics 
+
+This directory serves to aggregate all files necessary for multiplicity 
+and centrality calibration going from pp to Pb-Pb. It offers functionality 
+to perform simple slicing in percentiles, such as what is done in 
+proton-proton collisions, as well as the Glauber + analytical NBD fitter 
+used to perform anchoring and hadronic event distribution estimates in 
+nucleus-nucleus collisions.  
+
+## Files
+* `README.md` this readme
+* `CMakeLists.txt` definition of source files that need compiling
+* `multCalibrator.cxx/h` a class to do percentile slicing of a given histogram. Used for all systems. 
+* `multMCCalibrator.cxx/h` a class to perform data-to-mc matching of average Nch. 
+* `multGlauberNBDFitter.cxx/h` a class to do glauber fits. 
+* `multModule.h` a class to perform calculations of multiplicity and centrality tables for analysis. Meant to be used inside the main core service wagon 'multcenttable'. 
+

Common/Tools/Multiplicity/macros/README.md

@@ -0,0 +1,69 @@
+# Example multiplicity calibration macros
+
+You will find some example macros in this directory that will allow for the calculation of a glauber fit and the corresponding calibration histograms. 
+
+A simplified description of the procedure necessary to generate a 
+Glauber + NBD fit to a certain histogram is described below. 
+
+## Performing a Glauber + NBD fit
+
+### First step: calculation of Glauber MC sample 
+
+The Glauber + NBD model assumes that the multiplicity / amplitude
+distribution that one intends to fit can be described as coming 
+from a certain number N_{ancestor} of particle-emitting sources called 'ancestors'. Each ancestor is assumed to produce particles 
+according to a negative binominal distribution. Further, 
+the number of ancestors is assumed to be related to the basic 
+glauber quantities N_{part} and N_{coll} via the formula: 
+
+N_{ancestor} = f * N_{part}  + (1-f) N_{coll} 
+
+Usually, the value f is fixed at 0.8, and thus the range of
+N_{ancestors} is typically 0-700 in Pb-Pb collisions. 
+
+In order to allow for Glauber + NBD fitting, the correlation of 
+(N_{part}, N_{coll}) needs to be known. For that purpose, 
+a tree of Glauber MC needs to be generated using TGlauberMC using the relevant nuclei, which also involves choosing an appropriate nuclear profile. 
+
+Once TGlauberMC has been used to produce a tree with N_{part} and 
+N_{coll} values, their correlation needs to be saved to a 2D histogram that serves as input to the Glauber + NBD fitter used in 
+O2/O2Physics. This is done using the macro called `saveCorrelation.C` contained in this directory, which produces a file named `baseHistos.root`. The file `saveCorrelation.C` serves as 
+an example for the Pb-Pb case and minor adaptation may be needed
+in case other nuclei are to be used. 
+
+### Second step: execution of Glauber + NBD fitter 
+
+The fitting procedure is done via the macro ``. Because 
+the numerical fitting utilizes a convolution of a N_{ancestor} 
+distribution and NBD distributions, it will not be as fast 
+as typical one-dimensional fits: typical times of 10-100 seconds
+are not unusual. The macro will produce an output file 
+that contains: 
+
+* the original fitted histogram 
+* the actual glauber fit distribution, plotted all the way 
+to zero multiplicity 
+
+This output can then be used to extract percentiles. 
+
+### Third step: calculation of percentile boundaries 
+
+Once both the data and glauber fit distributions are known, 
+the next step is to create a 'stitched' data/glauber distribution in
+which the distribution at low multiplicities follows 
+the glauber fit and the distribution at higher multiplicities 
+follows the data. The multiplicity value in which the switch from 
+glauber to data takes place is called 'anchor point'. 
+
+Because of the fact that this 'stitching' procedure may need to be tuned and the actual glauber fit is slow, the stitching is done 
+in a separate macro called `runCalibration.C`. It is provided 
+in this directory as well and it is the third and last step in calculating percentile boundaries. The macro will printout some
+key boundaries as well as save an output file with the calibration. 
+
+*Bonus*: at the stage in which the glauber fit has been done 
+and all information is available, a de-convolution process
+can be used to calculate the average N_{part} and N_{coll} 
+in centrality slices. This functionality is also provided 
+in O2Physics as part of the `multGlauberNBDFitter` and the
+`runCalibration.C` macro can optionally also perform that 
+deconvolution. *Warning*: this procedure might take a mooment. 

Common/Tools/Multiplicity/macros/runCalibration.C

@@ -0,0 +1,150 @@
+void runCalibration( TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false){
+  TFile *file = new TFile(lInputFileName.Data(), "READ");
+  file->ls();
+  
+  TH1F *hData = (TH1F*) file->Get("hV0MUltraFine");
+  TH1F *hGlauberParameters = (TH1F*) file->Get("hGlauberParameters");
+  TH1F *hGlauberFitRange = (TH1F*) file->Get("hGlauberFitRange");
+  hData->SetName("hData");
+  TH1F *hStitched = (TH1F*) hData->Clone("hStitched");
+  TH1F *hFit = (TH1F*) file->Get("hGlauber");
+  
+  TCanvas *c1 = new TCanvas("c1", "", 800, 600);
+  c1->SetLeftMargin(0.17);
+  c1->SetBottomMargin(0.17);
+  c1->SetRightMargin(0.15);
+  c1->SetTopMargin(0.05);
+  c1->SetTicks(1,1);
+  c1->SetLogz();
+  c1->SetFrameFillStyle(0);
+  c1->SetFillStyle(0);
+  
+  
+  cout<<"Data bin width: "<<hData->GetBinWidth(1)<<endl;
+  cout<<"Fit bin width: "<<hFit->GetBinWidth(1)<<endl;
+  cout<<"Match range to use: "<<matchRange<<endl;
+  
+  //____________________________________________
+  double anchorPointFraction = anchorPointPercentage/100.f;
+  double anchorPoint = -1; // the anchor point value in raw
+  
+  //____________________________________________
+  // doing partial integration up to certain point for finding anchor point bin
+  for(int ii=1; ii<hData->GetNbinsX()+1; ii++){
+    // renormalize data curve
+    int bin1 = ii+1;
+    int bin2 = hData->FindBin( hData->GetBinLowEdge(ii+1) + matchRange + 1e-3 );
+    double matchRangeData = hData -> Integral( bin1, bin2);
+    double matchRangeFit = hFit -> Integral( bin1, bin2);
+    
+    // rescale fit to match in the vicinity of the region we're at
+    hFit->Scale(matchRangeData/matchRangeFit);
+    
+    double integralFit = hFit->Integral(1,ii);
+    double integralData = hData->Integral(ii+1,hData->GetNbinsX()+1);
+    double integralAll = integralFit+integralData;
+    
+    cout<<"at bin #"<<ii<<", integrated up to "<<hData->GetBinLowEdge(ii+1)<<" fraction above this value is: "<<integralData/integralAll<<endl;
+    anchorPoint = hData->GetBinLowEdge(ii+1);
+    
+    if(integralData/integralAll<anchorPointFraction) break;
+  }
+  
+  //____________________________________________
+  for(int ii=1; ii<hData->GetNbinsX()+1; ii++){
+    // renormalize data curve
+    if( hData->GetBinCenter(ii) < anchorPoint) hStitched->SetBinContent(ii, hFit->GetBinContent(ii));
+  }
+  
+  
+  cout<<"Anchor point determined to be: "<<anchorPoint<<endl;
+  cout<<"Preparing stitched histogram ... "<<endl;
+  
+  hFit->SetLineColor(kRed);
+  hStitched->SetLineColor(kBlue);
+  
+  hData->GetYaxis()->SetTitleSize(0.055);
+  hData->GetXaxis()->SetTitleSize(0.055);
+  hData->GetYaxis()->SetLabelSize(0.04);
+  hData->GetXaxis()->SetLabelSize(0.04);
+  hData->SetTitle("");
+  hData->Draw("hist");
+  hFit->Draw("hist same");
+  hStitched->Draw("hist same");
+  
+  //All fine, let's try the calibrator
+  multCalibrator *lCalib = new multCalibrator("lCalib");
+  lCalib->SetAnchorPointPercentage(100.0f);
+  lCalib->SetAnchorPointRaw(-1e-6);
+  
+  //Set standard Pb-Pb boundaries
+  lCalib->SetStandardOnePercentBoundaries();
+  
+  TString calibFileName = lInputFileName.Data();
+  calibFileName.ReplaceAll("glauberNBD", "calibration");
+  TFile *fileCalib = new TFile(calibFileName.Data(), "RECREATE");
+  
+  TH1F *hCalib = lCalib->GetCalibrationHistogram(hStitched, "hCalib");
+  
+  TCanvas *c2 = new TCanvas("c2", "", 800, 600);
+  c2->SetLeftMargin(0.17);
+  c2->SetBottomMargin(0.17);
+  c2->SetRightMargin(0.15);
+  c2->SetTopMargin(0.05);
+  c2->SetTicks(1,1);
+  //  c2->SetLogz();
+  c2->SetFrameFillStyle(0);
+  c2->SetFillStyle(0);
+  
+  hCalib->GetYaxis()->SetTitleSize(0.055);
+  hCalib->GetXaxis()->SetTitleSize(0.055);
+  hCalib->GetYaxis()->SetLabelSize(0.04);
+  hCalib->GetXaxis()->SetLabelSize(0.04);
+  hCalib->SetTitle("");
+  hCalib->Draw();
+  
+  
+  fileCalib->cd();
+
+  hData->Write();
+  hCalib->Write();
+  hStitched->Write();
+  hFit->Write();
+  
+  if(doNpartNcoll){
+    cout<<"Will now attempt to calculate % -> Np, Nc map..."<<endl;
+    
+    TProfile *hProfileNpart = new TProfile("hProfileNpart", "", 100, 0, 100);
+    TProfile *hProfileNcoll = new TProfile("hProfileNcoll", "", 100, 0, 100);
+    TH2F *h2dNpart = new TH2F("h2dNpart", "", 100, 0, 100, 500, -0.5f, 499.5f);
+    TH2F *h2dNcoll = new TH2F("h2dNcoll", "", 100, 0, 100, 3000, -0.5f, 2999.5);
+    
+    // Replay
+    multGlauberNBDFitter *g = new multGlauberNBDFitter("lglau");
+    TF1 *fitfunc = g->GetGlauberNBD();
+    
+    //Step 1: open the (Npart, Ncoll) pair information, provide
+    TFile *fbasefile = new TFile("basehistos.root","READ");
+    TH2D *hNpNc = (TH2D*) fbasefile->Get("hNpNc");
+    g->SetNpartNcollCorrelation(hNpNc);
+    g->InitializeNpNc();
+    
+    fitfunc->SetParameter(0, hGlauberParameters->GetBinContent(1));
+    fitfunc->SetParameter(1, hGlauberParameters->GetBinContent(2));
+    fitfunc->SetParameter(2, hGlauberParameters->GetBinContent(3));
+    fitfunc->SetParameter(3, hGlauberParameters->GetBinContent(4));
+    fitfunc->SetParameter(4, hGlauberParameters->GetBinContent(5));
+    
+    Double_t lMax = hData->GetBinLowEdge( hData->GetNbinsX() + 1);
+    
+    // uncomment if Np Nc needed -> Warning, slow!
+    g->CalculateAvNpNc( hProfileNpart, hProfileNcoll, h2dNpart, h2dNcoll, hCalib, 0, lMax );
+    
+    hProfileNpart->Write();
+    hProfileNcoll->Write();
+    h2dNpart->Write();
+    h2dNcoll->Write();
+  }
+
+  fileCalib->Write();
+}