]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/AliFragmentationFunctionCorrections.cxx
o First Version of TRDnSigma implementation (Xianguo) o still requires some catching...
[u/mrichter/AliRoot.git] / PWGJE / AliFragmentationFunctionCorrections.cxx
index 1770d3aa58dd7e3bf793e6f2e061ac1dc8ce058f..7d3a577ad3bf7f2b4a148dee17a608e1ab99e914 100644 (file)
 #include "AliFragmentationFunctionCorrections.h"
 #include <iostream> // OB TEST!!!
 #include <fstream>
+using std::cout;
+using std::endl;
+using std::cerr;
+using std::fstream;
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
@@ -72,7 +76,6 @@ AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
    ,fCorrBgr(0)
    ,fNCorrectionLevelsSinglePt(0)
    ,fCorrSinglePt(0)
-   ,fh1FFXiShift(0)
    ,fh1EffSinglePt(0)
    ,fh1EffPt(0) 
    ,fh1EffZ(0) 
@@ -142,7 +145,6 @@ AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const
   ,fCorrBgr(copy.fCorrBgr)    
   ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
   ,fCorrSinglePt(copy.fCorrSinglePt)   
-  ,fh1FFXiShift(copy.fh1FFXiShift)                  
   ,fh1EffSinglePt(copy.fh1EffSinglePt)
   ,fh1EffPt(copy.fh1EffPt)                      
   ,fh1EffZ(copy.fh1EffZ)                       
@@ -217,7 +219,6 @@ AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operat
     fCorrBgr                        = o.fCorrBgr;                      
     fNCorrectionLevelsSinglePt      = o.fNCorrectionLevelsSinglePt;
     fCorrSinglePt                   = o.fCorrSinglePt;
-    fh1FFXiShift                    = o.fh1FFXiShift;                  
     fh1EffSinglePt                  = o.fh1EffSinglePt;
     fh1EffPt                        = o.fh1EffPt;                      
     fh1EffZ                         = o.fh1EffZ;                       
@@ -275,8 +276,6 @@ AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
   if(fNJets)       delete fNJets;
   if(fNJetsBgr)    delete fNJetsBgr;
 
-  DeleteHistoArray(fh1FFXiShift);
-
   DeleteHistoArray(fh1EffPt);
   DeleteHistoArray(fh1EffXi);
   DeleteHistoArray(fh1EffZ );
@@ -746,8 +745,6 @@ void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const In
   fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
   AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
 
-  fh1FFXiShift = BookHistoArray();
-
   // eff histos
 
   fh1EffPt = BookHistoArray();
@@ -910,7 +907,7 @@ void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, c
 }
 
 //_____________________________________________________________________________________________________________________________________
-TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t jetPtSlice,  const Int_t type)
+TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t& jetPtSlice, const Int_t& type)
 { 
   // set histo bins for jet pt slice
   // if binning undefined for any slice, original binning will be used
@@ -1512,8 +1509,6 @@ void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString s
     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries())       fCorrFF[c]->GetZ(i)->Write();
     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries())      fCorrFF[c]->GetXi(i)->Write();
 
-    if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write();
-
     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries())       fCorrBgr[c]->GetZ(i)->Write();
     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries())      fCorrBgr[c]->GetXi(i)->Write();
@@ -2110,139 +2105,6 @@ void AliFragmentationFunctionCorrections::EffCorrBgr()
   }
 }
 
-//______________________________________________________________________
-void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
-{
-  // re-evaluate jet energy after FF corrections from dN/dpt distribution
-  // apply correction (shift) to dN/dxi distribution: xi = ln (pt/E) -> xi' = ln (pt/E') = ln (pt/E x E/E') = xi + ln E/E'
-  // argument corrlevel: which level of correction to be corrected/shifted to 
-
-  if(corrLevel>=fNCorrectionLevels){ 
-    Printf(" calc xi shift: corrLevel exceeded - do nothing");
-    return;
-  }
-
-  Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
-  Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
-  Double_t* deltaXi     = new Double_t[fNJetPtSlices];
-
-  for(Int_t i=0; i<fNJetPtSlices; i++){
-    
-    TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
-    TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
-
-    Double_t ptUncorr = 0;
-    Double_t ptCorr   = 0;
-
-    for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
-
-      Double_t cont   = histPtRaw->GetBinContent(bin);
-      Double_t width  = histPtRaw->GetBinWidth(bin);
-      Double_t meanPt = histPtRaw->GetBinCenter(bin);
-
-      ptUncorr += meanPt*cont*width;
-    }
-    
-    for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
-      
-      Double_t cont   = histPt->GetBinContent(bin);
-      Double_t width  = histPt->GetBinWidth(bin);
-      Double_t meanPt = histPt->GetBinCenter(bin);
-      
-      ptCorr += meanPt*cont*width;
-    }
-
-    jetPtUncorr[i] = ptUncorr; 
-    jetPtCorr[i]   = ptCorr;   
-  }
-
-  // ---
-
-  for(Int_t i=0; i<fNJetPtSlices; i++){
-
-    std::cout<<" jet pt bin "<<i<<" from "<<fJetPtSlices->At(i)<<" to "<<fJetPtSlices->At(i+1)
-            <<" jetPtUncorr "<<jetPtUncorr[i]<<" corrLevel "<<corrLevel<<" jet pt corr "<<jetPtCorr[i]<<" ratio "
-            <<jetPtCorr[i]/jetPtUncorr[i]<<std::endl;
-  } 
-
-  return; // for TEST!!!
-
-  // calc dXi from dN/dpt distribution : 
-  // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt 
-
-  for(Int_t i=0; i<fNJetPtSlices; i++){
-
-    Float_t jetPtLoLim = fJetPtSlices->At(i);
-    Float_t jetPtUpLim = fJetPtSlices->At(i+1);
-
-    Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
-    
-    Double_t ptUncorr = jetPtUncorr[i]; 
-    Double_t ptCorr   = jetPtCorr[i]; 
-
-    Double_t dXi = TMath::Log(ptCorr/ptUncorr);
-    
-    Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
-          ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
-    
-    deltaXi[i] = dXi;
-  }
-  
-  // book & fill new dN/dxi histos
-
-  for(Int_t i=0; i<fNJetPtSlices; i++){
-
-    TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
-            
-    Double_t dXi = deltaXi[i];
-
-    Int_t nBins  = histXi->GetNbinsX();
-    const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
-    Float_t binsVecNew[nBins+1];
-    
-    TString strName = histXi->GetName(); 
-    strName.Append("_shift");
-    TString strTit  = histXi->GetTitle();  
-
-    TH1F* hXiCorr; 
-
-    // create shifted histo ...
-
-    if(binsVec){ // binsVec only neq NULL if histo was rebinned before
-
-      for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;    
-      hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
-    }
-    else{ // uniform bin size
-      
-      Double_t xMin  = histXi->GetXaxis()->GetXmin();
-      Double_t xMax  = histXi->GetXaxis()->GetXmax();
-
-      xMin += dXi;
-      xMax += dXi;
-      
-      hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
-    }
-
-    // ... and fill
-
-    for(Int_t bin=1; bin<nBins+1; bin++){
-      Double_t cont = histXi->GetBinContent(bin);
-      Double_t err  = histXi->GetBinError(bin);
-      
-      hXiCorr->SetBinContent(bin,cont);
-      hXiCorr->SetBinError(bin,err);
-    }
-    
-    new(fh1FFXiShift[i]) TH1F(*hXiCorr);
-    delete hXiCorr;
-  }
-
-  delete[] jetPtUncorr;
-  delete[] jetPtCorr;
-  delete[] deltaXi;
-}
-
 //_____________________________________________________
 void AliFragmentationFunctionCorrections::SubtractBgr(Double_t sysErr)
 {
@@ -3009,21 +2871,21 @@ void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TSt
 
 //________________________________________________________________________________________________________________
 void AliFragmentationFunctionCorrections::WriteBgrJetSecCorr(TString strInfile, TString strBgrID, TString strID, TString strOutfile, 
-                                                            Bool_t updateOutfile, TString strOutDir)
+                                                            Bool_t updateOutfile, TString strOutDir,Double_t scaleFacBgrRec)
 { 
   // read task ouput from MC and write secondary correction - standard dir/list 
      
   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
   TString strlist = "fracfunc_" + strID;
     
-  WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID);
+  WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID,scaleFacBgrRec);
 }
 
 
 //___________________________________________________________________________________________________________________________________
 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
                                                          TString strOutfile, Bool_t updateOutfile, TString strOutDir,
-                                                         Bool_t writeBgr, TString strBgrID)
+                                                         Bool_t writeBgr, TString strBgrID,Double_t scaleFacBgrRec)
 {
   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
   // argument strlist optional - read from directory strdir if not specified
@@ -3298,9 +3160,12 @@ void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TSt
     hdNdzMCSecRecSsc[i]->SetNameTitle(strNameFFZSecRecSsc,"");
     hdNdxiMCSecRecSsc[i]->SetNameTitle(strNameFFXiSecRecSsc,"");
     
-    // normalize
+    // normalize 
     Double_t nJetsBin = fh1FFJetPtRecEffRec->Integral(binLo,binUp);
 
+    // scale fac for perp2 bgr
+    if(writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBin /= scaleFacBgrRec;
+
     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
@@ -4331,8 +4196,8 @@ void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile
   
 
   Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
-  Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
-  Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
+  Double_t binMinResponse[4] = { static_cast<Double_t>(loLimRecPt), static_cast<Double_t>(loLimTrackPt), static_cast<Double_t>(loLimGenPt), static_cast<Double_t>(loLimTrackPt)};
+  Double_t binMaxResponse[4] = { static_cast<Double_t>(upLimRecPt), static_cast<Double_t>(upLimTrackPt), static_cast<Double_t>(upLimGenPt), static_cast<Double_t>(upLimTrackPt)};
 
   const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
   
@@ -4735,7 +4600,7 @@ void AliFragmentationFunctionCorrections::dNdz2dNdxi()
 //________________________________________________________________________________________________________________
 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strIDGen,  TString strIDRec,  
                                                            TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
-                                                           TString strOutDir)
+                                                           TString strOutDir,Double_t scaleFacBgrRec)
 { 
   TString strdirGen  = "PWGJE_FragmentationFunction_" + strIDGen;
   TString strlistGen = "fracfunc_" + strIDGen;
@@ -4743,13 +4608,13 @@ void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, T
   TString strdirRec  = "PWGJE_FragmentationFunction_" + strIDRec;
   TString strlistRec = "fracfunc_" + strIDRec;
   
-  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kFALSE, "");
+  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kFALSE, "",scaleFacBgrRec);
 }
 
 //________________________________________________________________________________________________________________
 void AliFragmentationFunctionCorrections::WriteBgrBinShiftCorr(TString strInfile, TString strBgrID, TString strIDGen,  TString strIDRec,  
                                                               TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
-                                                              TString strOutDir)
+                                                              TString strOutDir,Double_t scaleFacBgrRec)
 { 
   TString strdirGen  = "PWGJE_FragmentationFunction_" + strIDGen;
   TString strlistGen = "fracfunc_" + strIDGen;
@@ -4757,14 +4622,14 @@ void AliFragmentationFunctionCorrections::WriteBgrBinShiftCorr(TString strInfile
   TString strdirRec  = "PWGJE_FragmentationFunction_" + strIDRec;
   TString strlistRec = "fracfunc_" + strIDRec;
   
-  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kTRUE, strBgrID);
+  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kTRUE, strBgrID,scaleFacBgrRec);
 }
 
 //___________________________________________________________________________________________________________________________________
 void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strdirGen, TString strlistGen, 
                                                            TString strdirRec, TString strlistRec, 
                                                            TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
-                                                           TString strOutDir,  Bool_t writeBgr, TString strBgrID)
+                                                           TString strOutDir,  Bool_t writeBgr, TString strBgrID, Double_t scaleFacBgrRec)
 {
  
   if((writeBgr && strBgrID.Length() == 0) || (!writeBgr && strBgrID.Length()>0) ){
@@ -5012,6 +4877,9 @@ void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, T
     Double_t nJetsBinGen = fh1FFJetPtMCGen->Integral(binLo,binUp);
     Double_t nJetsBinRec = fh1FFJetPtMCRec->Integral(binLo,binUp);
 
+    // scale fac for perp2 bgr
+     if(useRecPrim && writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBinRec /= scaleFacBgrRec;
+
     NormalizeTH1(hdNdptTracksMCGen[i],nJetsBinGen); 
     NormalizeTH1(hdNdzMCGen[i],nJetsBinGen); 
     NormalizeTH1(hdNdxiMCGen[i],nJetsBinGen); 
@@ -5791,7 +5659,8 @@ void AliFragmentationFunctionCorrections::ReadBinShiftCorrSinglePt(TString strfi
   }
  
 
-  if(fNHistoBinsPt) hBbBCorrPt = (TH1F*) hBbBCorrPt->Rebin(fNHistoBinsSinglePt,strNameBbBCorrPt+"_rebin",fHistoBinsSinglePt->GetArray());
+  if(fNHistoBinsPt && hBbBCorrPt) 
+    hBbBCorrPt = (TH1F*) hBbBCorrPt->Rebin(fNHistoBinsSinglePt,strNameBbBCorrPt+"_rebin",fHistoBinsSinglePt->GetArray());
 
   if(hBbBCorrPt) hBbBCorrPt->SetDirectory(0);