PWGLF/FORWARD/analysis2/AliForwardMultiplicityDistribution.cxx

   1 /**
   2  * @file   AliForwardMultiplicityDistribution.cxx
   3  * @author Christian Holm Christensen <cholm@master.hehi.nbi.dk>
   4  * @date   Thu Feb  7 01:03:52 2013
   5  *
   6  * @brief
   7  *
   8  * @ingroup pwglf_forward_multdist
   9  *
  10  */
  11 #include <TH1D.h>
  12 #include "AliForwardMultiplicityDistribution.h"
  13 #include "AliAODForwardMult.h"
  14 #include "AliAODCentralMult.h"
  15 #include "AliAODEvent.h"
  16 #include "TString.h"
  17
  18 ClassImp(AliForwardMultiplicityDistribution)
  19 #if 0
  20 ; // This is for Emacs - do not delete
  21 #endif
  22 //______________________________________________________________________
  23 AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution()
  24   : AliBaseAODTask(),
  25     fBins(),       // eta bin list
  26     fNBins(-1)     // multiplicity axis' runs from 0 to fNbins
  27 {
  28   //
  29   // Default Constructor
  30   //
  31 }
  32
  33 //_____________________________________________________________________
  34 AliForwardMultiplicityDistribution::AliForwardMultiplicityDistribution(const char* name)
  35   : AliBaseAODTask(name),
  36     fBins(),       // eta bin list
  37     fNBins(-1)    // multiplicity axis' runs from 0 to fNbins
  38 {
  39   //
  40   // Constructor
  41   //
  42 }
  43 //_____________________________________________________________________
  44 Bool_t AliForwardMultiplicityDistribution::Book()
  45 {
  46   TH1D* dndetaSumForward = new TH1D("dndetaSumForward","dndetaSumForward",
  47                                     200,-4,6);
  48   TH1D* dndetaSumCentral = new TH1D("dndetaSumCentral","dndetaSumCentral",
  49                                     200,-4,6);
  50   fSums->Add(dndetaSumForward);
  51   fSums->Add(dndetaSumCentral);
  52
  53
  54   TH1D* dndetaEventForward = new TH1D();
  55   TH1D* dndetaEventCentral = new TH1D();
  56   fSums->Add(dndetaEventForward);
  57   fSums->Add(dndetaEventCentral);
  58
  59   //Loop over all individual eta bins, and define their hisograms
  60   TIter next(&fBins);
  61   Bin * bin = 0;
  62   while ((bin = static_cast<Bin*>(next()))) {
  63     bin->CreateOutputObjectss(fSums, fNBins);
  64   }
  65   return true;
  66 }
  67
  68
  69 //_____________________________________________________________________
  70 Bool_t AliForwardMultiplicityDistribution::Event(AliAODEvent& aod)
  71 {
  72   //
  73   // User Exec
  74   //
  75   AliAODForwardMult* AODforward = GetForward(aod);
  76   if (!AODforward) return false;
  77
  78   AliAODCentralMult* AODcentral = GetCentral(aod);
  79   if (!AODcentral) return false;
  80
  81   TH2D& forward = AODforward->GetHistogram();
  82   TH2D& central = AODcentral->GetHistogram();
  83
  84   TH1D* dndetaSumForward   = (TH1D*)fSums->FindObject("dndetaSumForward");
  85   TH1D* dndetaSumCentral   = (TH1D*)fSums->FindObject("dndetaSumCentral");
  86   TH1D* dndetaEventForward = 0;//(TH1D*)fList->FindObject("dndetaEventForward");
  87   TH1D* dndetaEventCentral = 0;//(TH1D*)fList->FindObject("dndetaEventCentral");
  88   TH1D* normEventForward   = 0;
  89   TH1D* normEventCentral   = 0;
  90
  91   // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
  92   dndetaEventForward = forward.ProjectionX("dndetaForward",1,forward.GetNbinsY(),"");
  93   dndetaEventCentral = central.ProjectionX("dndetaCentral",1,central.GetNbinsY(),"");
  94   dndetaSumForward->Add(dndetaEventForward);
  95   dndetaSumCentral->Add(dndetaEventCentral);
  96
  97   // underflow eta bin of forward/central carry information on whether
  98   // there is acceptance. 1= acceptance, 0= no acceptance
  99   //
 100   // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
 101   normEventForward   = forward.ProjectionX("dndetaEventForwardNorm",0,0,"");
 102   normEventCentral   = central.ProjectionX("dndetaEventCentralNorm",0,0,"");
 103
 104   Double_t VtxZ = AODforward->GetIpZ();
 105
 106
 107   // loop over all eta bins, and fill multiplicity histos
 108   TIter next(&fBins);
 109   Bin * bin = 0;
 110   while ((bin = static_cast<Bin*>(next()))) {
 111     bin->Process(dndetaEventForward, dndetaEventCentral,
 112                  normEventForward,   normEventCentral, VtxZ);
 113   }
 114
 115   // ACHTUNG ACHTUNG ACHTUNG! Serious memory leak here!
 116   // We should do
 117   //
 118   // delete dndetaEventForward;
 119   // delete dndetaEventCentral;
 120   // delete normEventForward;
 121   // delete normEventCentral;
 122   //
 123   // and remove them from the output list - to prevent the memory leak
 124
 125   return true;
 126
 127 }
 128
 129 //______________________________________________________________________
 130 #if 1
 131 const Char_t*
 132 #else
 133 TString
 134 #endif
 135 AliForwardMultiplicityDistribution::FormBinName(Double_t etaLow, Double_t etaHigh)
 136 {
 137   //
 138   //  Form name of eta bin
 139   //
 140   TString sLow(TString::Format("%+03d",int(10*etaLow)));
 141   TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
 142   sLow.ReplaceAll("+", "plus");
 143   sLow.ReplaceAll("-", "minus");
 144   sHigh.ReplaceAll("+", "plus");
 145   sHigh.ReplaceAll("-", "minus");
 146  #if 0
 147   TString* combined= new TString();
 148   combined->Append(TString::Format("%s_%s", sLow.Data(), sHigh.Data()));
 149   return combined->Data();
 150 #else
 151   static TString tmp;
 152   tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
 153   return tmp.Data();
 154 #endif
 155 }
 156
 157 //=====================================================================
 158 AliForwardMultiplicityDistribution::Bin::Bin()
 159   : TNamed(),
 160     fEtaLow(-1111),     // low eta limit
 161     fEtaHigh(-1111),    // high eta limit
 162     fHist(0),           // multiplicity distribution hist
 163     fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
 164     fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
 165     fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
 166     fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
 167     fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
 168     fHistMinus10(0),    // multiplicity distribution hist scaled down with 10%
 169     fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
 170     fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
 171     fAcceptance(0),     // histogram showing the 'holes' in acceptance.
 172     fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 173 {
 174   //
 175   // Default Constructor
 176   //
 177 }
 178 //_____________________________________________________________________
 179 AliForwardMultiplicityDistribution::Bin::Bin(Double_t etaLow, Double_t etaHigh)
 180   : TNamed("", ""),
 181     fEtaLow(etaLow),    // low eta limit
 182     fEtaHigh(etaHigh),  // high eta limit
 183     fHist(0),           // multiplicity distribution hist
 184     fHistPlus05(0),     // multiplicity distribution hist scaled up with 5%
 185     fHistPlus075(0),    // multiplicity distribution hist scaled up with 7.5%
 186     fHistPlus10(0),     // multiplicity distribution hist scaled up with 10%
 187     fHistMinus05(0),    // multiplicity distribution hist scaled down with 5%
 188     fHistMinus075(0),   // multiplicity distribution hist scaled down with 7.5%
 189     fHistMinus10(0),    // multiplicity distribution hist scaled down with 10%
 190     fHistPlusSys(0),    // multiplicity distribution hist scaled up with the event uncertainty
 191     fHistMinusSys(0),   // multiplicity distribution hist scaled down with the event uncertainty
 192     fAcceptance(0),     // histogram showing the 'holes' in acceptance.
 193     fVtxZvsNdataBins(0) // VtxZ vs. number of data acceptance bins (normalised to the eta range)
 194 {
 195   //
 196   // Constructor
 197   //
 198   const char* name = AliForwardMultiplicityDistribution::FormBinName(fEtaLow,fEtaHigh);
 199
 200   SetName(name);
 201   SetTitle(TString::Format("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
 202
 203 }
 204 //_____________________________________________________________________
 205 void AliForwardMultiplicityDistribution::Bin::CreateOutputObjectss(TList* cont,  Int_t max)
 206 {
 207   //
 208   // Define eta bin output histograms
 209   //
 210   TList* out = new TList;
 211   out->SetName(GetName());
 212   cont->Add(out);
 213
 214   fHist            = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
 215   fHistPlus05      = new TH1D("multPlus05", GetTitle(), max, -0.5, max-.5);
 216   fHistPlus075     = new TH1D("multPlus075", GetTitle(), max, -0.5, max-.5);
 217   fHistPlus10      = new TH1D("multPlus10", GetTitle(), max, -0.5, max-.5);
 218   fHistMinus05     = new TH1D("multMinus05", GetTitle(), max, -0.5, max-.5);
 219   fHistMinus075    = new TH1D("multMinus075", GetTitle(), max, -0.5, max-.5);
 220   fHistMinus10     = new TH1D("multMinus10", GetTitle(), max, -0.5, max-.5);
 221   fHistPlusSys     = new TH1D("multPlusSys", GetTitle(), max, -0.5, max-.5);
 222   fHistMinusSys    = new TH1D("multMinusSys", GetTitle(), max, -0.5, max-.5);
 223   fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins", "VtxZ vs dataAcceptance/etaRange;z-vtz;dataAcceptance/etaRange", 20, -10,10, 130,0,1.3);
 224   fAcceptance      = new TH2D("Acceptance","Acceptance;#eta;z-vtx",200,-4, 6 , 20,-10,10);
 225
 226   out->Add(fHist);
 227   out->Add(fHistPlus05);
 228   out->Add(fHistPlus075);
 229   out->Add(fHistPlus10);
 230   out->Add(fHistMinus05);
 231   out->Add(fHistMinus075);
 232   out->Add(fHistMinus10);
 233   out->Add(fHistPlusSys);
 234   out->Add(fHistMinusSys);
 235   out->Add(fVtxZvsNdataBins);
 236   out->Add(fAcceptance);
 237
 238 }
 239
 240 //_____________________________________________________________________
 241 void AliForwardMultiplicityDistribution::Bin::Process(TH1D* dndetaForward, TH1D* dndetaCentral,
 242                                                       TH1D* normForward,   TH1D* normCentral, Double_t VtxZ)
 243 {
 244   //
 245   // Process single eta bin
 246   //
 247
 248   // Diagnostics on event acceptance
 249   Int_t    first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
 250   Int_t    last  = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
 251   Double_t acceptanceBins=0;
 252   for(Int_t n= first;n<=last;n++){
 253     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
 254       acceptanceBins++;
 255     }
 256     fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ), 1);
 257     if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
 258       fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(VtxZ),10);
 259   }
 260
 261   //std::cout << "% of bins covered in eta-range: " << acceptanceBins/(last-first+1) << std::endl;
 262   fVtxZvsNdataBins->Fill(VtxZ, (Double_t)acceptanceBins/(last-first+1));
 263
 264
 265   Double_t c     = 0;
 266   Double_t e2    = 0;
 267   for (Int_t i = first; i <= last; i++) {
 268     Double_t cForward  = 0;
 269     Double_t cCentral  = 0;
 270     Double_t e2Forward = 0;
 271     Double_t e2Central = 0;
 272     if (normForward->GetBinContent(i) > 0) {
 273       cForward  = dndetaForward->GetBinContent(i);
 274       e2Forward = dndetaForward->GetBinError(i) * dndetaForward->GetBinError(i);
 275     }
 276     if (normCentral->GetBinContent(i) > 0) {
 277       cCentral  = dndetaCentral->GetBinContent(i);
 278       e2Central = dndetaCentral->GetBinError(i) * dndetaCentral->GetBinError(i);
 279     }
 280     Double_t cc  = 0;
 281     Double_t ee2 = 0;
 282     // if overlap between central/forward, use average of the two
 283     if (cCentral > 0 && cForward > 0) {
 284       cc  = 0.5 * (cForward  + cCentral);
 285       ee2 = 0.5 * (e2Forward + e2Central);
 286     }
 287     else if (cCentral > 0) {
 288       cc  = cCentral;
 289       ee2 = e2Central;
 290     }
 291     else if (cForward > 0) {
 292       cc  = cForward;
 293       ee2 = e2Forward;
 294     }
 295     c  += cc;
 296     e2 += ee2;
 297   }
 298
 299   //Systematic errors from here
 300
 301   Double_t fmd=0;
 302   Double_t spd=0;
 303   Double_t overlap=0;
 304
 305   // number of eta bins in fmd, spd and overlap respectively
 306   for(Int_t i = first;i<=last;i++){
 307     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
 308       fmd++;
 309     if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
 310       overlap++;
 311     if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
 312       spd++;
 313   }
 314
 315   Double_t sysErrorSquared = 0;
 316   // estimate of systematic uncertainty on the event multiplicity - hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or Casper Nygaard phd theses.
 317   Double_t fmdSysError= 0.08;
 318   Double_t spdSysError= 0.067;
 319
 320   Double_t total = 0;
 321   total= fmd+ spd + overlap;
 322   // Combined systematc event uncertainty, by weighting with the number of eta-bins of fmd-only, spd-only and the overlap.
 323   if(total){
 324     sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+ spd*TMath::Power(spdSysError,2)+
 325                       0.5*overlap*TMath::Power(fmdSysError,2)+ 0.5*overlap*TMath::Power(spdSysError,2))/total;
 326   }
 327
 328   //Strangeness correction, assumed to be 2.0 +- 1% (syst),  taken from HHD and CN thesis
 329   c=0.98*c;
 330
 331
 332   // correction for missing material description. Correction is based on separate PbPb analysis of difference between nominel and displaced vertices
 333   if(fEtaHigh<1.55&&fEtaHigh>1.45)
 334     c=0.98*c;
 335   if(fEtaHigh<2.05&&fEtaHigh>1.95)
 336     c=0.963*c;
 337   if(fEtaHigh<2.45&&fEtaHigh>2.35)
 338     c=0.956*c;
 339   if(fEtaHigh>2.95)
 340     c=0.955*c;
 341
 342   fHist->Fill(c);
 343   fHistPlusSys->Fill((1+TMath::Sqrt(sysErrorSquared))*c);
 344   fHistMinusSys->Fill((1-TMath::Sqrt(sysErrorSquared))*c);
 345   fHistPlus05->Fill(1.05*c);
 346   fHistPlus075->Fill(1.075*c);
 347   fHistPlus10->Fill(1.1*c);
 348   fHistMinus05->Fill(0.95*c);
 349   fHistMinus075->Fill(0.925*c);
 350   fHistMinus10->Fill(0.9*c);
 351
 352
 353 }
 354
 355
 356
 357
 358 //_____________________________________________________________________
 359 //
 360 //
 361 // EOF